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ABSTRACT 

We exploit the deep and extended far-infrared data-sets (at 70, 100 and 160 /im) of 
the Herschel GTO PACS Evolutionary Probe (PEP) Survey, in combination with the 
HERschel Muhi-tiered Extragalactic Survey (HerMES) data at 250, 350 and 500 ^m, 
to derive the evolution of the rest-frame 35- //m, 60-^m, 90-/xm, and total infrared (IR) 
luminosity functions (LFs) up to z~4. We detect very strong luminosity evolution 
for the total IR LF (Ljr cx(l-hz)^-^^=^"-^° up to z~2, and cx(l-|-z)i-62±o.5i 2<z<4) 
combined with a density evolution (cx(l-|-z)~°'^^**''^^ up to z^l and oc(l+z)^'^'^^*''-^'* 
at l<z^4). In agreement with previous findings, the IR luminosity density (pir) 
increases steeply to z^l, then flattens between z~l and z~3 to decrease at z^3. 
Galaxies with different SEDs, masses and sSFRs evolve in very different ways and 
this large and deep statistical sample is the first one allowing us to separately study 
the different evolutionary behaviours of the individual IR populations contributing to 
piR. Galaxies occupying the well established SFR~stellar mass main sequence (MS) are 
found to dominate both the total IR LF and pm at all redshifts, with the contribution 
from off-MS sources (>0.6 dex above MS) being nearly constant (~20% of the total 
Pir) and showing no significant signs of increase with increasing z over the whole 
0.8<2<2.2 range. Sources with mass in the range lO^log(M/M0)^ll are found to 
dominate the total IR LF, with more massive galaxies prevailing at the bright end of 
the high-z {^2} LF. A two-fold evolutionary scheme for IR galaxies is envisaged: on the 
one hand, a starburst-dominated phase in which the SMBH grows and is obscured by 
dust , is followed by an AGN-dominated phase, then evolving toward a local elliptical. 
On the other hand, moderately star-forming galaxies containing a low-luminosity AGN 
have various properties suggesting they are good candidates for systems in a transition 
phase preceding the formation of steady spiral galaxies. 

Key words: cosmology: observations - galaxies: active ~ galaxies: evolution - galax- 
ies: luminosity function - galaxies: starburst - infrared: galaxies. 



1 INTRODUCTION 

Understanding the origin and growth of the galaxies we ob- 
serve today is one of the main problems of current cosmol- 
ogy. The luminosity function (LF) provides one of the fun- 
damental tools to probe the distribution of galaxies over 
cosmological time, since it allows us to assess the statisti- 
cal nature of galaxy formation and evolution. When com- 
puted at different redshifts, the LF constitutes the most di- 
rect method for exploring the evolution of a galaxy popula- 
tion, describing the relative number of sources of different 
luminosities counted in representative volumes of the Uni- 
verse. The LF computed for different samples of galaxies 
can provide a crucial comparison between the distribution 
of different galaxy types, i.e. galaxies at different redshifts, 
in different enviroments or selected at different wavelengths. 

It has now become clear that we cannot understand 
galaxy evolution without accounting for the energy absorbed 



by dust and re-emitted at longer wavelengths (e.g, Genzel 
I&: Cesarsky 20001, in the infrared (IR) or sub-millimetre 
(sub-mm). Dust is responsible for obscuring the ultraviolet 
(UV) and optical light from galaxies: since star-formation 
occurs within dusty molecular clouds, far-IR and sub-mm 
data, where the absorbed radiation is re-emitted, are es- 
sential for providing a complete picture of the history of 



* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and 
with important participation from NASA 
f E-mail; carlotta.gruppioniaoabo.inaf.it 



star-formation through cosmic time, which is one of the fun- 
damental instruments we have to reconstruct how galaxies 
have evolved since their formation epoch. For these reasons, 
extragalactic surveys in the rest-frame IR represent a key 
tool for understanding galaxy formation and evolution. 

Surveys of dust emission performed with the former 
satellites exploring the Universe in the mid- and far-IR do- 
main, i.e. the Infrared Astronomical Satellite {IRAS;\Neuge-\ 



bauer 19841 and the Infrared Space Observatory {ISO; 
Kessler et al. 1996 1, allowed the first studies of the IR-galaxy 



LF at z<0.3 (Saunders et al. 19901 and z<l (Pozzi et al 



20041, respectively. With Spitzer 24- /im data, it was pos- 
sible to study the evolution of the mid-IR LF up to z~2 
(e. g.|Le Floc'h et al.|200'5l [Caputi et al.|2007[ [Rodighiero et] 



al]|2010a[), while, even with the deepest Spitzer Space Tele- 
scope ( [Werner et al]|2004[ ) 70 -fim data, only z~l-1.2 could 
be reached in t he far-IR ([Magnelli et al.|p5o9l |Patel et ah] 
2012[) ) - though pagneUi et al.|(|2011[) reached z~2 through 



stacking. Since the rest-frame IR spectral energy distribu- 
tions (SEDs) of star-forming galaxies and AGN peak at 60- 
200 fim, to measure their bolometric luminosity and evolu- 
tion with z we need to observe in the far-IR/sub-mm regime. 
However, the detection of large numbers of high-z sources at 
the peak of their IR SED was not achievable before the Her- 
schel Space Observatory ( Pilbratt et al.|2010 l, due to source 
confusion and/or low detector sensitivity, and our knowledge 
of the far-IR luminosity function in the distant Universe is 
still affected by substantial uncertainties. Ground-based and 
balloon-borne observations in the mm/sub-mm range, prob- 
ing the evolution of the most distant (2 J; 2) and luminous 
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dusty galaxies, have so far been limited to the identification 
of sources at the very bright end of the luminosity function 
(e.g., Chapman et al.] '2005). All of these works detected 
strong evolution in both luminosity and/or density, indi- 
cating that IR galaxies were more luminous and/or more 
numerous in the past. Strong observational evidence of high 
rates of evolution for IR galaxies has been obtained also 
through the detection of a large amount of energy contained 



in the Cosmic Infrared Background (CIRB; Hauser & Dwek 



2001 1, and the source counts from several deep cosmological 



surveys (from 15 /xm to 850 ^m) largely exceeding the no- 
evolution expectations (e.g. Small et al. |1997[ [Elbaz et al. 



T999| JPapovich et al.|2004| [Bethermin et al.|2010||Marsden 



et al.|20lT l. Both the CIRB and the source counts require a 
strong increase in the IR energy density between the present 
time and 2~l-2. At higher redshifts the total emissivity 
of IR galaxies is poorly constrained, due to the scarcity of 
Spitzer galaxies at z>2, the large spectral extrapolations to 
derive the total IR luminosity from the mid-IR (see e.g. 



El- 



|baz et al.|20T0l|Nordon et al.|2010| and |Nordon et al.|2012 



for 

descriptions of the failure, at least at z>1.5, of previous to- 
tal IR luminosity extrapolations from the mid-IR, although 
we must note that this failure mainly affects luminosity- 
dependent methods like, e.g., that of Chary &: ElbazpOOl I 
and the incomplete information on the ^-distribution of sub- 
mm sources (Chapman et al. 2005). 

Herschel, with its 3.5-m mirror, is the first telescope 
which allows us to detect the far-IR population to high red- 
shifts (2:~4-5) and to derive its rate of evolution through 
a detailed LP analysis. The new extragalactic surveys pro- 
vided by Herschel in the far-IR/sub-mm domain, like the 
wide and shallow Herschel- ATLAS ( Eales et al.|2010 Dunne 



et al. 20111, the complementary Herschel Multi-tiered Ex- 
tragalactic Survey (HerME S; |01iver et alT]|2012| ) and PACS 



Evolutionary Probe (PEP; |Lutz et al. 2011 1 covering the 
most popular cosmological fields, and the deep, pencil beam. 



Herschel-GOOHS project ( jElbaz et al.||2011[ ), will be cru- 
cial to assess galaxy and AGN evolution in the IR at z>2. 
They will give us the opportunity to study in detail the 
population of IR galaxies and their evolution with cosmic 
time since the Universe was about a billion years old. In 
particular, the Photodetector Array Camera & Spectrome- 



ter {PACS; Poglitsch et al. 20101, with its high sensitivity 



and resolution at 70-/im, WO-fim and 160-^m, is the best 
suited instrument to detect faint IR sources by overcoming 
the source confusion and blending problems that affected 
the previous far-IR missions. 

This is the first of two papers aiming at deriving the 
far- and total IR LFs from the Herschel PACS-\- Spectral 



and Photometric Imaging Receiver {SPIRE; jGriffin et al. 



20101 data obtained within the PEP and HerMES extra- 



galactic survey projects. In the present paper, we derive the 
rest-frame 60-/im, 90-/im and total IR (8-1000 /im) 

LFs from a sample selected at PACS 70, 100 and 160 /im 
wavelengths in the GOODS (GOODS-S and GOODS-N), 
Extended Chandra Deep Field South (ECDFS) and COS- 
MOS areas. We use the full 70-500 PACS-f SPIRE data 
to determine Lir and SED properties of the PACS selected 
sources. In a related paper, Vaccari et al. (in prep.) derive 
rest-frame 100-, 160- and 250-/.im and total IR LFs for a 
SPIRE selected sample. In addition, a third work aimed at 
studying the total IR LF based on the 24-/im selected sam- 



ple, using all the PEP-^HerMES data in the COSMOS field, 
is ongoing (Le Floc'h et al., in preparation). 

PEP is one of the major Herschel Guaranteed Time ex- 
tragalactic key-projects, designed specifically to determine 
the cosmic evolution of dusty star-formation and of the IR 
luminosity function. It is structured as a "wedding cake", 
based on four different layers covering different areas to dif- 
ferent depths at 100 and 160 fim (in the GOODS-S field also 
at 70 /im), from the large and shallow COSMOS field to the 
deep, pencil beam GOODS-S field. PEP includes the most 
popular and widely studied extragalactic fields with exten- 
sive multi- wavelength coverage available, in particular deep 
optical, near-IR and Spitzer imaging and spectroscopic and 
photometric redshifts: COSMOS; Lockman Hole; Extended 
Grot h Streep (EGS); ECDFS; GOODS-N ; an d GOODS-S 
(see iBerta et al.||2010| |Berta et al.||2011| and |Lutz et al 



[2011 



for a detailed description of the fields and observa- 
tions). Coordinated observations of the PEP fields at 250, 
350 and 500 fim with SPIRE have been obtained by the Her- 



MES Survey (Oliver et al. 2012 1. HerMES, analogously to 



PEP but extending to a much wider area, is a legacy pro- 
gramme designed to map a set of nested fields (~380 deg^ 
in total) of different sizes and depths, using SPIRE (at 250, 
350 and 500 /.im), and PACS (at 100 and 160 /^m, shallower 
than PEP), with the widest component of 270 deg'^ with 
SPIRE alone. In the fields covered by PEP, the two surveys 
are closely coordinated to provide an optimized sampling 
over wavelength. 

In [Gruppioni et al.] ( |2010[ ) we started to determine the 
evolution with redshift of the galaxy and AGN LF in the far- 
IR domain by exploiting the PEP data obtained in GOODS- 
N by the PEP Science Demonstration Programme (SDP). 
Here we extend the analysis to the wider and shallower 
fields - COSMOS and ECDFS - and to the deepest field - 
GOODS-S - observed by PEP, and we also take advantage 
of the HerMES sub-mm data in the same fields to derive 
improved SED classifications and accurate total IR lumi- 
nosities for our sources. This allows us to have statistically 
significant samples of IR galaxies at different redshifts and 
over a broad range of luminosities, to make a detailed study 
of the LF at several z intervals, all the way from z—Q to 2;~4. 
The measure of the total IR luminosity obtained by inte- 
grating the SEDs, well constrained over the entire mid- and 
far-IR domain (and also in the sub-mm thanks to the avail- 
able SPIRE data), allows us to derive the total IR LF and 
its evolution directly from far-IR data for unbiased samples 
selected at wavelengths close to the peak of dust emission. 
Moreover, the availability of deep multi-wavelength cata- 
logues in the PEP fields is crucial for analysing the SEDs, 
obtaining k-corrections and total IR luminosities, and classi- 
fying the PEP sources into different IR populations, in order 
to separately study their LFs and evolutionary behaviour. 
This is the first study ever based on such a statistically wide 
and deep far-IR sample, to be able to provide LFs for differ- 
ent IR populations of galaxies and AGN. Here the evolution 
of the far- and total IR LFs (and luminosity density, here- 
after pir) are derived up to unprecedented high redshifts 
(~4) both globally (e.g. for all the populations together) 
and separately, for each SED class. 

Despite the abundance of information available in the 



literature about the stellar mass function (MF; Fontana et 



|al.|2004||Pozzetti et al.|2010|[llbert et al.|20l"0l|Domingue^ 
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Sanchez et al. 20111, very little is known about the corre- 



sponding total IR LFs and star-formation rate (SFR) den- 
sities at different masses (an attempt based on Spitzer data 
was made by Perez-Gonzalez et al. 2005J . From stellar MF 
studies one finds a clear increase with z of the relative frac- 
tion of massive (log(M/M0>ll) star-forming objects, start- 
ing to contribute significantly to the massive-end of the MFs 
at z>l ( [Fontana et al.||2004[ [llbert et al.|2010[ ). Their evo- 
lution and contribution to the total SFR history is however 
still uncertain, since only few studies have tried to recon- 
struct the evolution of the SF history of massive objects 
from optical/near-IR or mid-IR surveys ( Juneau et al.|2005{ 
Pcrcz-G onzalez et al.|[2005l [Santini et al.||2009| [Fontanot et| 
al.^^2012J but none from far-IR selected surveys (providing 
a more direct indicator of the galaxy SF activity). In this 
work, we have derived the IR luminosity function and den- 
sity in three different mass ranges (from log(M/MQ)=8.5 
to log(M/MQ)=12), extending previous studies (limited to 
2=1.8-2 for the most massive galaxies) to z~4. 
Finally, our PEP data-sets have allowed us to quantify the 
relative contribution of the two main modes of star forma- 
tion (a relatively steady one in disk-like galaxies, defining 
a tight SFR-stellar mass sequence, and a starburst mode in 
outliers) to the total IR LF and in three redshift inter- 
vals (0.8<2<1.25, 1.25<2<1.8 and 1.8<2<2.2) and to test 
the SED-classes belonging to each mode. 

The paper is structured as follows. The PEP Survey 
with the far-IR and multi-wavelength data, together with 
the SED characterisation and redshift distribution of the 
PEP sources, are described in Sect. 2. The LFs (rest-frame 
35-/im, 60-^m, 90-fim and total IR), their evolution (derived 
for different SED-classes, mass and specific star-formation 
rate intervals) are discussed in Sect. 3. In Sect. 4 we present 
the number and IR luminosity densities of the difi'erent 
galaxy types, while in Sect. 5 we discuss our results. In 
Sect. 6 we present our conclusions. 

Throughout this paper, we use a Chabrier initial mass 
function (IMF) and we assume a ACDM cosmology with 
Ho = 71 km s"^ Mpc"\ = 0.27, and SIa = 0.73. 



2 THE DATA 

The PEP fields where we computed the LFs are: COSMOS, 
2 deg^ observed down to 3ct depths of ~5 mjy and 10.2 
mjy at 100 /im and 160 /xm, respectively; ECDFS, ~700 
arcmin^ down to 3a depths of ~4.5 mJy and 8.5 mJy at 
100 fim and 160 fim, respectively; GOODS-N, ~300 arcmin^ 
to 3 and 5.7 mJy at 100 pim and 160 /^m, respectively; and 
GOODS-S, ~300 arcmin^ to 1.2 mJy, 1.2 mJy and 2.4 
mJy at 70 /im, 100 fim and 160 fim, respectively. Our ref- 
erence samples are the blind catalogues at 70 (in GOODS-S 
only), 100 and 160 /xm to the 3a level, which contain 373 
(all in GOODS-S), 7176 (GOODS-S: 717, GOODS-N: 291, 
ECDFS: 813, COSMOS: 5355) and 7376 (GOODS-S: 867, 
GOODS-N: 316, ECDFS: 688, COSMOS: 510 5) sources at 
70, 100 and 160 /xm, respectively. We refer to |Berta et al.| 
( |2010[ ) and [Berta et aL\ ( |2011[ ) for a detailed description of 
the data catalogues and source counts. 



2.1 Multi-wavelength Identification 

The PEP fields benefit from an extensive multi-wavelength 
coverage. We have therefore associated our sources to the an- 
cillary catalogues by means of a multi-band likelihood ratio 
technique (Suth erland fc Saunders^ 1992 : Cilicgi eral.|2001 1, 
starting from the longest available wavelength (160 /xm, 
PACS) and progressively matching 100 /xm (PACS), 70 /xm 
(PACS, GOODS-S only) and 24 /xm {Spitzer /MIPS). In the 
GOODS-S field, we have associated to our PEP sources 
the 24:-fim catalogue by Magnelli et al. (20091, that we 



have matched with the optical-l-near-IR-l-IRAC MUSIC cat- 



alogue of Grazian et al. (20061, revised by Santini et al, 
(2009]), which includes spectroscopic and photometric red- 
shifts. To maximise the fraction of identifications, we lim- 
ited our study to the area covered by the MUSIC cata- 
logue (~196 arcmin^), obtaining 233, 468 and 492 sources at 
70, 100 and 160 fim, respectively, with fiux density greater 
than the flux limits reported above (all with either spec- 
troscopic or photometric redshifts). In the GOODS-N field. 



as described in Berta et al. (20101, Berta et al. (20111 and 



Gruppiorii_ et al. (20101, a PSF-matched multi- wavelength 
catalogu^was created, including photometry from the far- 
UV [GALEX) to the mid-IR [Spitzer). As in GOODS-S, 
to maximise the identifications, we limited our study in 
GOODS-N to the area covered by the ACS (~150 arcmin^), 
obtaining 176 and 186 sources with fiux density greater than 
the flux limit at 100 and 160 /xm, respectively (all with red- 
shifts). We have matched our sources in the ECDFS with the 
multi-wavelength Su rvey by Yale-Chile (MUSYC) by |Car- 
damone et al. (20101, obtaining 687 sources at 100 /im and 
625 sources at 160 /im (578 and 547 with redshifts, ~45% 
spectroscopic). Finally, in COSMOS, we have matched our 
catalogue with the deep 24-/xm sample of |Le Floc'h et ah 
( 2009 1 and with the IRAC-based catalogue of llbert et al. 



(20101, including optical and near-IR photometry and pho- 



tometric redshifts. After the removal of PEP sources within 
flagged areas of the optical and/or IRAC COSMOS cata- 
logues, we ended up with two catalogues consisting of 4110 
and 4118 sources, with flux densities Js5.0 and ^10.2 mJy 
at 100 and 160 /xm respectively (3817 and 3849 with either 
spectroscopic or photometric redshifts). Throughout this pa- 
per and speciflcally for the SED fits described in Section [2.2[ 
we adopt these spectroscopic or rest-frame UV to near-IR 
photometric redshifts for the various fields. 



The HerMES extragalactic survey ( Oliver et al. 2012 1 



performed coordinated observations with SPIRE at 250, 350 
and 500 /im in the same fields covered by PEP. In particular, 
in HerMES a prior source extraction was performed using 



the method presented in Roseboom et al. (2011), based on 



MIPS-24 /tm positions. The 24-/im sources used as priors for 
SPIRE source extraction are the same as those associated 
with our PEP sources through the likelihood ratio technique. 
We have therefore associated the HerMES sources with the 
PEP sources by means of the 24- /im sources matched to both 
samples. For most of our PEP sources (~87 per cent) we 
found a >3a SPIRE counterpart in the HerMES catalogues. 



^ publicly available at 

http: //www.mpe .mpg. de/ir/Research/PEP/public_data_releases . 
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Figure 1. Observed rest-frame SEDs of the PEP sources (black dots) divided by population (as shown in the plot) and normalized to 
the Xg-band. The more representative templates for each SED-class have been overplotted in different colours. 



2.2 Galaxy Classification 

We made use of all the available multi-wavelength data to 
derive tlie SEDs of our PEP sources, whicli we interpreted 
and cl assified by pe rforming a x'^ fit (using tiie Le Phare 



codtj^ [Arnouts et~aLj(2062] and [llbert et aL][2006[) wit h tiie 
semi-empirical template library of PoUetta et al. ( 2007 1, rep- 



resentative of different classes of IR galaxies and AGN. To 
this library we added some templates modified in their far- 
IR part to better reproduce the observed Herschel data (see 
Gruppioni et al.]|2010| ), and three starburst templates from 
Rieke et al. (2009). If required to improve the fit, different 



^ available at 

http : //www. cf ht .hawail . edu/~arnouts/LEPHARE/lephare .html 



extinction values (i?(B_v) from 0.0 to 0.5) have been ap- 
plied to the templates, by letting the code free to choose 
the most suitable extinction curve. The considered set of 
templates included SEDs of elliptical galaxies of different 
ages, lenticular, spirals (from Sa to Sdm), starburst galax- 
ies (SB), type 1 QSOs, type 2 QSOs, Seyferts, LINERs and 
composite ULIRGs (containing both starburst and obscured 
AGN component), in the wavelength range between 0.1 and 
1000 /im. The latter templates, are empirical ones created to 
reproduce the SEDs of the heavily obscured AGN. Two of 
these SEDs (the broad absorption-line QSO Markarian 231 
(Bcrta 20 05|) and the S eyfert 2 galaxy IRAS 19254-7245 
South (Berta et al.|2003| )) are similar in shape, containing a 
powerful starburst component, mainly responsible for their 
far-IR emission, and an AGN component that contributes to 
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Table 1. SED Classification of the PEP Sources 



nciQ 


spiral 


St arburst 






R\jri 1 


or Ryjii yiutij 


017— APM/'Q-i-,-; T-al 
or JlLilM V,opii a.J. J 


GOODS-S 


N 


N 


N 


N 


N 


N 


N 


70 fim 


53 


22 


142 






26 


116 




(23%) 


(9%) 


(61%) 


(2%) 


(5%) 






100 /xm 


117 


60 


250 


10 


31 


54 


96 




(25%) 


(13%) 


(53%) 


(2%) 


(7%) 






160 fim 


123 


55 


277 


11 


26 


73 


204 




(25%) 


(11%) 


(56%) 


(2%) 


(6%) 






GOODS-N 


N 


N 


N 


N 


N 


N 


N 


100 /^m 


68 


20 


78 


7 


3 


21 


57 




(39%) 


(11%) 


(44%) 


(4%) 


(2%) 






160 fim 


67 


21 


85 


10 


3 


21 


64 




(36%) 


(11%) 


(46%) 


(5%) 


(2%) 






ECDFS 


N 


N 


N 


N 


N 


N 


N 


100 /^m 


253 


49 


245 


8 


23 


83 


162 




(44%) 


(9%) 


(42%) 


(1%) 


(4%) 






160 fim 


233 


49 


231 


12 


22 


99 


132 




(43%) 


(9%) 


(42%) 


(2%) 


(4%) 






COSMOS 


N 


N 


N 


N 


N 


N 


N 


100 Atm 


1637 


232 


1689 


76 


183 


580 


1109 




(43%) 


(6%) 


(44%) 


(2%) 


(5%) 






160 fim 


1483 


243 


1847 


103 


173 


777 


1070 




(39%) 


(6%) 


(48%) 


(3%) 


(4%) 






TOTAL 


N 


N 


N 


N 


N 


N 


N 


100 fim 


2075 


361 


2262 


101 


240 


738 


1424 




(41%) 


(7%) 


(45%) 


(2%) 


(5%) 






160 fim 


1906 


368 


2440 


136 


224 


970 


1470 




(38%) 


(7%) 


(48%) 


(3%) 


(4%) 







- and dominates - the mid-IR (Farrah et al. 2003), reproduc- 
ing the SEDs of "obscured" AGN regardless of their optical 



spectra (i.e. broad or narrow lines in the optical; Gruppioni 
|et al.|2008 1. Hereafter, we will refer to this class of templates 



and to the sources reproduced by them as to type 2 AGN 
(AGN2). Three other empirical templates, reproducing the 
observed SEDs of nearby ULIRGs containing an obscured 
AGN (i.e. IRAS 20551-4250; IRAS 22491-1808; NGC 6240) 
have been associated to the Seyfert 1.8/2, LINER ones, since 
they all contain an AGN, but this AGN does not dominate 
the observed energetic output at any wavelength (from UV 
to far-IR/sub-mm), showing up just in the range where the 
host galaxy SED has a minimum (i.e. the mid-IR). The AGN 
in these objects is either obscured or of low luminosity. We 
refer to this class as to star-forming galaxies containing an 
AGN (SF-AGN), since their IR luminosity is largely domi- 
nated by star-formation. 

In our analysis, we make the basic assumption that the SED 
shapes seen at low redshifts are also able to represent the 
higher redshift objects. In any case, to further increase the 
range of SEDs in the fit, we have applied additional extinc- 
tion with different extinction curves to our templates. All 
SED fits adopt fixed spectroscopic or photometric redshifts 
described in Section |2T] 

The template library used to fit our data contains a finite 
number of SEDs (38), representative of given classes of lo- 



cal infrared objects, which do not vary with continuity from 
one class to another (there are large gaps in the param- 
eter space). Therefore, the quality of the fit depends not 
only on the photometric errors, but also on the template 
SED uncertainties. For this reason, in our fitting procedure, 
in addition to the photometric errors on data, we need to 
take into account also the uncertainties due to the tem- 
plate SEDs discretisation and additional extinction. To do 
this, we have proceeded as described in detail by |Gruppi-] 



oni et al. (20081 and summarised as follows. First, we have 



run Le Phare on our PEP SEDs considering the nominal 
errors from catalogues, computing the distributions of the 

(•S'objcct ~ S'tcmplatG)band/(o-objoct)band ValuCS iu each of the 

considered photometric band (where Sobjoct and cTobject are 
the fiux density and the relative error of the source, and 
•S'tcmpiatc the flux dcusity of the template in the considered 
band), iteratively increasing the photometric errors until we 
have obtained a Gaussian distribution with (t~1. This cor- 
responds to reduced distributions peaked around 1 (as 
expected in the case of good flt). With the new photomet- 
ric uncertainties (on average, significantly increased mainly 
in the optical/near-IR and SPIRE bands), we have run Le 
Phare on our sources for the second time, obtaining what 
we have taken as the final SED-fitting results. 

The majority of our PEP sources are reproduced by 
templates of normal spiral galaxies (spiral), SB galax- 
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ies (starburst), and Seyfert2/1.8/LINERS/ULIRGs+AGN 
(SF-AGN), although different classes prevail at different red- 
shifts and luminosities. The spiral SEDs show no clear signs 
of enhanced SF or nuclear activity (see Fig.[l]|, the far-lR 
bump being characterised by relatively cold dust (rdust~20 
K). On the other hand, SB templates are characterised by 
warmer (Tdust~40-45 K), more pronounced far-IR bumps 
and significant UV extinction, indicative of intense star- 
formation activity. Templates of star-forming galaxies con- 
taining either a low-luminosity or obscured AGN (SF-AGN) 
are characterised by a "flattening" in the 3-10 /xm spectrum 
(suggesting detection of an AGN in the wavelength range 
where the host galaxy SED has a minimum) and a far-IR 
bump dominated by star-formation, which is intermediate 
(in terms of both energy and Tdust) between spirals and SBs. 
Although they can be considered as star-forming galaxies at 
the wavelengths relevant to this work, we prefer to refer to 
them as SF-AGN throughout the paper, to keep in mind that 
they probably contain an AGN, whose presence, though not 
dominant in the far-IR, might be very important for analy- 
sis in other bands (e.g., in the X-rays or the mid-IR). 
We note that the well-studied high-z (~2.3) SED of the 
strongly lensed sub-mm gala:xy SMM J2135-0102, known as 
"the Cosmic Eyelash" (Ivison et al. 2010 Swinbank et al 



20101, best-fitted with our procedure by an extincted (E(B- 
V)~0.2) IRAS 22491-1808 template (though rather poorly 
in the near-IR), does not represent the bulk of our popu- 
lation at high-2 (>1.5), whose SEDs are indeed well repro- 
duced by our library of templates. 

In fact, the considered template set provides very good 
fits to the SEDs of our PEP sources. In Fig. [l] we show the 
rest-frame SEDs (black dots) of the PEP sources belonging 
to the different "broad" SED classes (spiral, starburst, 
SF-AGN, AGNl and AGN2), compared to the template SEDs of 
those classes normalised to the A"s-band flux density. In Ta- 
ble[l]we report the fraction of sources belonging to each SED 
class: we find that in all the fields ~41(38) per cent of the 
100(160)-/im sources are reproduced by a spiral template 
SED, 7(7) per cent with a starburst template SED, 45(48) 
per cent with a SF-AGN template SED, 2(3) per cent with an 
AGN2 SED and 5(4) per cent with an AGNl SED. We note that 
the fraction of SF-AGN derived in this work is in agreement 
with results from mid-IR spectroscopy (with Spitzer-IRS) 
of local star-forming galaxies from the SINGS sample by 
( |2007) , who found that ~50 per cent of local 
galaxies (though of lower luminosities than ours) do harbour 
low- luminosity AGN (of LINER or Seyfert types). Recently, 
Sajina et al. (20121 found an even higher fraction (~70 per 
cent) of objects hosting an AGN in the mid-IR {Spitzer 24- 
/im) selected samples (~23 per cent AGN-dominated and 
~47 per cent showing both AGN and starburst activity). 
However, since the far-IR SED of the SF-AGN is dominated 
by star-formation and at these wavelengths resembles either 
starburst or spiral galaxy templates, we have also divided 
the SF-AGN class into SF-AGN (SB) and SF-AGN (Spiral) sub- 
classes, based of their far-IR/near-IR colours (e.g. Sioo/Si.e) 
and SED resemblance (apart from the rest-frame mid-IR 
flattening, which is detected in all of the SF-AGN SEDs). 
Specifically, galaxies best-fitted by the Seyfert2/1.8 tem- 
plates (either the original ones from Polletta et al.|[2"007 • 
those modified by [Gruppioni et al ' 



2010[ ) have been clas- 



NGC 6240, IRAS 20551-4250 or IRAS 22491-1808 templates 
have been classified as SF-AGN (SB). The number of sources 
belonging to the former and the latter sub-classes are also 
reported in Table [l] as additional information. 

2.3 Redshift Distribution 

A large number of spectroscopic redshifts have been mea- 
sured in the GOODS, ECDFS and COSMOS regions. In 
the GOODS-S and ECDFS area a collection of more than 



5000 spectroscopic redshifts are available (( 


I!ristiani et al. 


20001 ICroom et al. |2001 Bunker et al. |2003 


|Dickinson et 


al.||2004| IStanway et al. |2004 |Strolger et al. 


2004| ISzokoly 


et al. 2004 van der Wei et al. 2004 


Doherty et al. 2005| 


Le Fevre et al. 2005| Mignoli et al. 


20051 Vanzella et al. 


2008; Popesso et al. 2009; Santini et al. 2009; Balestra et 



,al.^2010, Cooper et al._,2012) . In the GOODS-N area more 
than 2000 spectroscopic redshifts come from various obser 
vatio ns ([Cohen et al.||2000| |Wirth et aT]|2004[ |Cowie et aT 



2004i [Barger et al.||2008||. Finally, in COSMOS we could 



use a collection of ~3000 spectroscopic redshifts from ei- 
ther the public zCOSMOS bright database or the non-public 



al 



20091 



zCOSMOS deep database (Lilly, S.J. et al. 2007 Lilly et 



For the PEP sources without spectroscopic red- 
shift available, we have adopted the photometric redshifts 
derived from multi-wavelength (UV to near-IR) photome- 
try by different authors in the different fields, as mentioned 
in Section |2]T] In the GOODS-S field the MUSIC photo- 



metric redshift catalogue (Grazian et al. 2006 Santini et 



al. 2009 1 provided photo- zs for most of our PEP sources 



without spectroscopic data, while in the GOODS-N field. 



photo- as were obtained by Berta et al. (20101 for almost all 



the PEP sources within the ACS area. The ICardamone eti 
[al] ( [20T0| ) and [Ubert et al!] ( [2009| catalo gues provided pho- 
tometric redshifts for a large fraction of the PEP sources in 
the ECDFS and COSMOS areas, respectively. When con- 
sidering both the spectroscopic and photometric redshifts, 
in our PEP fields the redshift incompleteness is very low. In 
particular, in the GOODS-S field we have either a spec-z or 
a photo- 2: for ~100 per cent of the PEP sample within the 
MUSIC areas (~80 per cent spectroscopic, though most of 
them lie at z<2.5; see jBerta et al.|2011[ ). In the GOODS-N 
field we have a redshift completeness of ~100 per cent of 
sources (70 per cent spectroscopic) within the ACS area. In 
the ECDFS and COSMOS fields we have a redshift com- 
pleteness of 88 per cent and 93 per cent respectively (45 per 
cent and 40 per cent spectroscopic). 

The uncertainty in the photometric redshifts has been 
evaluated by means of a comparison with the available 
spec-zs by the different authors providing photo-a cata- 



logues in the PEP fields. In particular, Berta et al. (20111 
have compared the photometric and the available spectro- 
scopic redshifts in GOODS-S, GOODS-N and COSMOS, 
finding a fraction of outliers, defined as objects having 
Az/(l+2spcc)>0.2, of ~2 per cent for sources with a PACS 
detection. Most of these outliers are sources with few pho- 
tometric points available, or SEDs not well reproduced by 
the available templates. The median absolute deviation of 
the Az/{1+Zapec) distribution in the three fields analised 



sified as SF-AGN (Spiral) , while galaxies best-fitted by the 



by Berta et al. (20111 is 0.04 for the whole catalogue, and 
0.0 38 for PACS-detected objects. In GOODS-S, [Grazian et| 
al. ( 2006 1 found an excellent agreement between photomet- 
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Figure 2. Redshift distributions of 100-/xm (top) and 160-/^m (bottom) sources in the four PEP fields (first row from top, GOODS- 
S; second row, GOODS-N; third row, ECDFS; bottom row, COSMOS) to different Hmiting fluxes (as shown in the plot). The redshift 
distributions of the five different populations have been plotted in different colours (green, spiral; cyan, starburst; red, SF-AGN; magenta, 
AGN2; blue, AGNl) and compared to the total distribution (black solid histogram). The line-filled dashed histograms shown in the spiral 
and starburst panels represent the redshift distributions of the SF-AGN_Spiral and SF-AGN_SB sub-classes, respectively. 
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ric and spectroscopic redshifts over the fully accessible red- 
shift range 0<2:<6 {a[Az/{l + z)]~0.045), with a very lim- 
ited number of catastrophic errors. In COSMOS, [Ilbert et al.| 
(20101 estimated the photometric redshift uncertainties of 



Table 2. Average redshifts in the PEP Survey fields 



their catalogue matched with the COSMOS photo-z 

multi- wavelength catalogue of Ilbert et al. (20091, finding 
a[Az/{l + 2:)]=0.008 (and <1 per cent of catastrophic fail- 
ures) at i+0<22.5, cr[A2/(l-|-z)]=0.011 at 22.5<iXB<24 and 
a[Az/{l + z)]=0m3 at 24<jJ^B<25. In the ECDFS, by com- 
paring non-X-ray sources with high-quality spectroscopic 
redshifts, [Cardamone et al.|2010| found a [ Az/ ( 1 + z)] =0.008 
to z~1.2 =0.027 at 1.2^z^3.7 and =0.016 at z>3.7. Note 
that we have checked all the z>2.5 photometric redshifts 
through Le Phare, assigning the Le Phare derived value in 
case of significant disagreement with that from the cata- 
logue (though most resulted in very good agreement). The 
fractions of spectroscopic redshifts in the 2.5<2:<3.0 interval 
amount to just ~6% and ~25% in COSMOS and GOODSS 
respectively, dropping to ~4% and ~6% at 3.0<2:<4.2. We 
note that from our comparison between photo- and spec-z 
(when available) we find a general good agreement in all 
fields. 

Photometric redshift errors may, in principle, affect the 
shape of the luminosity function at the bright end: by scat- 
tering objects to higher redshifts they make the steep fall- 



off at high luminosities appear shallower (e.g. Drory et al. 
2003). To study the impact of redshift uncertainties on the 
inferred infrared LF, we have performed Monte Carlo sim- 
ulations, as discussed in detail in Section 13.31 It is indeed 



very difficult to estimate the effect of catastrophic failures 
at z>2.5, where mainly photometric redshifts are available 
and very little reliable spectroscopic data can be used to 
validate them. Moreover, for the limited high-z samples 
with spectroscopic information, different results are found 
in the different fields: i.e., in GOODS-S [Berta et al.| ( |2011 l 
found about 25 per cent of catastrophic failures for PACS 
detected sources above z~2, with the tendency to have a 
higher than real photometric redshift, while in COSMOS 
the catastrophic failures (20 per cent) found by [Ilbert et| 



al. (20091 for MIPS selected sources at 2<2;<3 were mostly 
for photo-2's smaller than spectroscopic ones. In addition to 
that, sometimes high-z spectroscopic redshifts can be even 
more uncertain than photometric ones and great care must 
be taken when selecting spec-zs for comparison (i.e. we need 
to choose those with high quality flags). For all these rea- 
sons, we limited our analysis of photo-z uncertainties to the 
Monte Carlo simulations described in Section |3.3[ without 
trying to derive uncertainties also due to catastrophic fail- 



ures. In Section 3.2 we also note that the different (far-IR) 
photo-z approach of Lapi et al. ( 20111) produces consistent 



LF results in the common part of parameter space. 

The median redshift of the 70-pim sample in GOODS-S 
is 2med(70)=0.67 (the mean is (2)70=0.86), while those of the 
100- and 160-/im samples are different in each of the fields, 
given the different flux density depths reached by PEP in 
each area. In Table [2] we report the median and the mean 
redshifts found for the different flelds (and for the combined 
sample) at the different selection wavelengths. As expected, 
GOODS-S reaches the highest redshifts, while the surveys in 
the COSMOS and ECDFS flelds are shallower and sample 
lower redshifts. On average, the 160-^m selection favours 



field 


70 fim 


100 Atm 


160^m 




^mcd 




^med 




^med 




GOODS-S 


0.67 


0.86 


0.85 


1.07 


0.98 


1.16 


GOODS-N 






0.73 


0.84 


0.84 


0.94 


ECDFS 






0.66 


0.76 


0.69 


0.85 


COSMOS 






0.59 


0.74 


0.70 


0.88 


ALL FIELDS 


0.67 


0.86 


0.64 


0.78 


0.73 


0.91 



higher redshifts than the 100-/im one (see also Berta et al. 
|20flp . 

In Fig. [2] we show the redshift distribution of the PEP 
sources selected at 100 /im and 160 /im in the four different 
flelds. The black solid histogram is the total redshift distri- 
bution in the fleld (one for each row of the plot), while the 
fllled histograms in different colours represent the redshift 
distributions of the different populations (green, spiral; 
cyan, starburst; red, SF-AGN; magenta, AGN2; blue, AGNl). 
The line-flUed dashed histograms shown in the spiral and 
starburst panels represent the redshift distributions of the 
SF-AGN(Spiral) and SF-AGN(SB) sub-classes, respectively. 
In addition to the principal redshift peak, in GOODS-S a 
secondary peak centred at 2~2 is clearly visible at both 
100 and 160 fim. A similar result has been shown and dis- 



cussed also by Berta et al. (20111, while an extensive anal- 
ysis of PACS GOODS-S large-scale structure at 2=2-3 and 
of a z~2.2 filamentary overdensity have been presented by 



Magliocchetti et al. (20111 



3 THE LUMINOSITY FUNCTION 

The sizes and depths of the PEP samples are such as to al- 
low a direct and accurate determination of the far-IR LF in 
several redshift bins, from 2~0 up to z~4. PEP+HerMES is 
the unique Herschel survey to allow such analysis over such a 
wide redshift and luminosity range, sampling both the faint 
and bright ends of the far- and total IR LFs with suflicient 
statistics. Because of the redshift range covered by PEP, 
we would need to make signiflcant extrapolations in wave- 
length when computing the rest-frame LFs at any chosen 
wavelengths. In order to apply the smallest extrapolations 
for the majority of our sources, we choose to derive the far- 
IR LFs at the rest-frame wavelengths corresponding to the 
median redshift of each sample. Given the median redshift 
of the 70- Aim sample in GOODS-S (~0.67, see Table [2}, we 
use that sample to derive the rest-frame luminosity function 
at 35 /im. With the 100- and 160-/im PEP samples (whose 
median redshifts are ~0.64 and 0.73 respectively), we de- 
rived the rest-frame LFs at 60 and 90 /^m. Note that, given 
the excellent multi-wavelength coverage available for most 
of our sources (thanks also to the HerMES data available 
in all the PEP flelds and providing reliable counterparts for 
most of our PEP sources), their SEDs are very well deter- 
mined from the UV to the sub-mm. The extrapolations are 
therefore well constrained by accurately defined SEDs, even 
at high redshifts (i.e. at 2~3.5 the rest-frame 90-/xm lumi- 
nosity corresponds to AobscrvGd~400 //m, which is still in the 
range covered by HerMES). 
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Figure 3. Rest-frame 35 /xm Luminosity Function estimated with the 1/Vmax method from the PEP 70-/jm sample in the GOODS-S 
field (purple open circles) and independently in the four PEP fields from the 100-fim selected samples (red filled circles, GOODS-S; green 
filled triangles, GOODS-N; orange filled stars, ECDFS; blue filled squares, COSMOS). The error-bars in the data points represent the 
Poissonian uncertainties. For comparison, we also plot the deter mination of|Magnelli et ah] | |2011[ ) at 1.2<z<2.5, shown as green open 
squares, and the double power-l aw fit of|Magnelli et al.| | |2009[ | and jMagnelli et al.| ( |2011[ l, shown as a green dashed line. The black dashed 
line is the z=0 determination of |Magnelli et al.| l |2009^ . 



3.1 Method 

The LFs are derived using the 1/Knax method (Schmidt 
1968). This method is non-parametric and does not require 
any assumptions on the LF shape, but derives the LP di- 
rectly from the data. We have first derived the LFs in each 
field separately, in order to check for consistency and to test 
the role of cosmic variance. Successively, we have made use 
of the whole data-sets to derive the monochromatic and to- 



tal IR LFs, by means of the Avni & Bahcall ( 1980( 1 method 
for coherent analysis of independent data-sets. We have di- 
vided the samples into different redshift bins, over the range 
^ z ^ 4, selected to be almost equally populated, at least up 
to z^2.5. In each redshift bin we have computed the comov- 
ing volume available to each source belonging to that bin, 
defined as 'Knax^Kniax — K„,i„, where jZmax is the minimum 
between the maximum redshift at which a source would still 
be included in the sample given the limiting flux of the sur- 
vey (different for each fleld) and the upper boundary of the 
considered redshift bin, while Zmin is just the lower bound- 
ary of the considered redshift bin. 

When combining the four samples, we have constructed 
a complete sample over the whole GOODS-S+GOODS- 
N-fECDFS(-GOODS-S)-HCOSMOS region, including aU 



the observed objects (see details in Section 3.2 1. The depth 
of the sample is not constant throughout the region, but 
an object with a given flux density (included in the list of 



observed objects irrespective of the identity of its parent 
sample) can a priori be found in one (or more) region if 
its redshift is ^S2:^'ix('S'iimit) of that region (e.g. sources de- 
tected in the COSMOS area are detectable over the whole 
joint area, while the fainter sources detected in GOODS-S 
are detectable in GOODS-S only). The maximum volume of 
space which is available to such an object to be included in 
the joint sample is then defined by 
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where V^^^^^ (with fld=GS, GN, E, C corresponding to 
GOODS-sT'gOODS-N, ECDFS and COSMOS, respec- 
tively) is the comoving volume available to each source in 
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Figure 4. Rest-frame 60 /im Luminosity Function estimated independently from the 100-^m selected samples in the four PEP fields 
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derived both at 0.0<2< 0.2 and at 0.2<2<0.4). The cyan open triangles show the PEP 60-Atm LF from the SDP data in GOODS-N by 
|Gruppioni et alT] | |2010[ l (the redshift bins are not exactly the same). The grey dashed line is the |Saunders et al.| l |19Sjo[ | local LF. 



that field, in a given redshift bin, while fifld is the solid angle 

subtended by that field sample on the sky. 

For each luminosity and redshift bin, the LF is given by: 



1 



AL ^ Wi X Vni 



(2) 



where Knax.i is the comoving volume over which the i-th 
galaxy could be observed, AL is the size of the luminosity 
bin, and Wi is the completeness correction factor of the i-th 
galaxy. These completeness correction factors are a combi- 
nation of the completeness corrections given by|Berta et al 



(20101 and Berta et al. (20111, derived as described in 



et al. 



Lutz 



(20111, to be applied to each source as function of 



its flux density, together with a correction for redshift in- 
completeness (for the ECDFS and COSMOS only, see Sec- 
tion 2.3 1. Since, as mentioned in 2.3 the redshift incomplete- 



ness in the COSMOS and ECDFS areas is independent on 
PACS flux density, in these fields we have applied the cor- 
rections regardless of the source luminosity and redshift (i.e. 
by muftiplying ^{L,z) by 1.07 and 1.14 in COSMOS and 
ECDFS, respectively). However, the redshift incompleteness 
does not affect our conclusions, since > 95 per cent of all our 
sources do have a redshift. 

Uncertainties in the infrared LF values depend on pho- 
tometric redshift uncertainties. To quantify the effects of 
the uncertainties in photometric redshifts on our luminosity 



functions, we performed a set of Monte Carlo simulations, 
as described in Section |3^ 



3.2 The Rest-Frame 35-, 60- and 90- /im 
Luminosity Function 

By following the method described above, we have derived 
the 35- /^m, 60- /im and 90- /im rest-frame LFs from the 70-fim 
(in GOODS-S only), 100-/im and 160-/im samples, respec- 
tively. In order to check the consistency between the cata- 
logues and the effects of cosmic variance, we have first de- 
rived the monochromatic LFs in each field separately. Note 
that, since the 70-fim data are available in the GOODS-S 
field only, to have a better sampling of the LF especially at 
the bright-end, we have also computed the rest-frame 
LF from the 100-pim samples in the four fields and compared 
them with that obtained from the 70- /im sample (see Fig|3] 
Table 3). The agreement between the two derivations is very 
good, implying correct extrapolations in wavelength due to 
the good and complete SED coverage. 

We have divided the samples into seven redshift bins: 
0.0<2s;0.4; 0.4<2^0.8; 0.8<ziil.2; 1.2<z^l.8; 1.8<2^2.5; 
2.5<2^3.5; and 3.5<z ^4.5. The results of the computa- 
tion of our 35-fim (reported in Table 3), 60-fim and 90-/im 
LFs are shown in Figs. [3] |4] and [5j respectively. The LFs in 
the four different fields appear to be consistent with each 
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Figure 5. Rest-frame 90 /im Luminosity Function estimated independently from the 160-^m selected samples in the four PEP fields 
with the 1/Vmax method (with the same symbols as in Fig. [4]l . Th e diagonal crosses represent the local LF of |Serjeant et ar| l |2004[ |, 
the brown dots with error bars are the local LF of jSedgwick et al.| | |2011| in the AKARI Deep Field, while the pink open circles are 
the 100-/im LF derivat ion of |Lapi et al] | |2011[ | in the H- ATLAS fields. The cyan open triangles show the PEP 90-fim LF computed by 
[Gruppioni et al~] l |2010[ | from the SDP data in the GOODSN (not exactly the same redshift bins). 



other within the error bars (±lo-) in most of the common 
luminosity bins. The COSMOS and GOODS-S Surveys are 
complementary, with the faint end of the LFs being mostly 
determined by data in GOODS-S, and the bright end by 
COSMOS data. After having checked the field-to-field con- 
sistency, we have combined the 100- and 160-/im samples in 
all fields, obtaining the global rest-frame 60- and 90-/im LFs 
(reported in Tables 4 and 5, respectively). 

The data from each field in each z-hin have been plotted 
(and considered in the combination) only in the luminosity 
bins where we expect our sample to be complete, given that 
at fainter luminosities not all galaxy types are observable 
(depending on their SEDs; Ilbert et al. 2004). For compar- 
ison, we overplot the LFs at 35 /im from |Magnelli et al.| 



( 2009[) and [Magnelh et al.| ( |2011[ l, the local LFs at 6 ^m 
from [Saunders et al.| ( 1990|) and thos e at 90 /xm from [Ser-[ 
jeant et al. (2004|) and Sedgwick et al. ( 2011 1 and at 100 fim 



from Lapi et al. (20111, respectively. The comparison be- 



tween the 35-fj.m LF, derived from the 70-/im PEP sample 



in GOODS-S, and the resuhs of Magnelli et al. ( 2009 1 and 



Magnelli et al. ( 2011 1, based on a 24- /im prior extraction and 



stacking analysis on Spitzer maps, shows very good agree- 
ment, both with the data and with the double power-law fit. 
The 1.8<2<2.5 PEP LF is consistent within ±la with the 
[Magnelli et al.| ( [2011[ ) data points, though the power-law fit 
at bright L35 (>10 Lq) seems to be slightly lower than our 



data. At z>2.5 no comparison data from Spitzer are avail- 
able, while our LF derivation can provide hints of evolution 
at the bright end of the LF. In the common redshift intervals 
(between 2~1 and 3.5), our 90-/im LF is in very good agree- 
ment with the WO-fj,m |Lapi et al. (20111 derivation from 
the H- ATLAS survey (although their redshift bins are some- 
what diflferent than ours: 1.2<2<1.6; 1.6<z<2.0; 2.0<z;<2.4; 
and 2.4<z<4.0) and with the previous PEP-SDP derivation 
( [Gruppioni et al.|2010 K The consistency between our 90- /xm 
LFs and the Lapi et al. (20111 ones (derived from a differ- 
ent sample, using a different template SED to fit the data 
and a rest-frame far-IR based method to derive photometric 
redshifts) gives us confidence that, at least up to z^3.5, our 
computation is not significantly affected by incompleteness 
or by photo-z uncertainties. 



3.3 The Total Infrared Luminosity Function 

We integrate the best-fit SED of each source over 
8^5 Arost^ 1000 /im to derive the total IR luminosities {Lm — 
L[8-1000/im]) in 11 redshift bins (0.0-0.3; 0.3-0.45; 0.45- 
0.6; 0.6-0.8; 0.8-1.0; 1.0-1.2; 1.2-1.7; 1.7-2.0; 2.0-2.5; 2.5- 
3.0; and 3.0-4.2), selected to be almost equally populated, 
at least up to 0~2.5. Our approach is similar to that of 
other studies based on mid-IR selected galaxy samples (e.g. 
[Le Floc'h eraL][2005l [Rodighiero et al.pOlOal [Magnelli et] 
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Figure 6. Total IR Luminosity Function estimated independently from the 160-fim samples in the four PEP fields (with the same 
symbols as in Fig.[4]l. The pink and sky-blue shaded areas represent the range of values derived with 20 iterations by allowing a change 
in photo-2 in the GOODS-S and COSMOS fields, respectively. 



al.|2011 1, but this is the first time that the SEDs have been 
accurately constrained by sufliciently deep data in the far- 
IR/sub-mm domain and not simply extrapolated from the 
mid-IR to the longer wavelengths or from average flux den- 
sity ratios. 

As mentioned in Section[23] we have studied the impact 
of redshift uncertainties on the total IR LF by performing 
Monte Carlo simulations. As test cases, we used the COS- 
MOS and GOODS-S samples (which are basically defining 
the bright and faint ends of the LF in an almost complemen- 
tary way), and we checked the effect on the total IR luminos- 
ity function. We iterated the computation of the total IR LF 
by each time varying the photometric redshift of each source 
(assigning a randomly selected value, according to a Gaus- 
sian distribution, within the 68 per cent confidence interval). 
Each time, we then recomputed the monochromatic and to- 
tal IR luminosities, as well as the Knax value, but we did 
not perform the SED-fitting again, keeping the previously 
found best-fit template for each object (the effect on the k- 
correction is not relevant in the far-IR wavelength range). 
The results of this Monte Carlo simulation are reported in 
Fig. [6j where we show the total IR LFs derived in each of 
the PEP fields independently: the red and blue filled cir- 



cles represent the estimates of the GOODS-S and COSMOS 
LFs, with their range of values derived with 20 iterations by 
allowing a change in photo-z represented by the pink and 
sky-blue shaded areas, respectively. The comparison shows 
that the effect of the uncertainty of the photometric red- 
shifts on the error bars is slightly larger than the simple 
Poissonian value and affects mainly the lower and 

the higher redshift bins (especially at low and high luminosi- 
ties). Using these Monte Carlo simulations, we find no ev- 
ident systematic offsets caused by the photometric redshift 
uncertainties (see Fig. [6|. This is due to the very accurate 
photometric redshifts available in these fields. For the total 
uncertainty in each luminosity bin in GOODS-S and COS- 
MOS, we have therefore assumed the dispersion given by the 
Monte Carlo simulations (as shown in Fig.[7|. We note that 
at the higher z (>2.5), where we must rely on a majority 
of photometric redshifts, the true uncertainties (taking into 
account also catastrophic failures or incompleteness effects) 
might be larger than derived with simulations. The unavail- 
ability of a "true" comparison sample (i.e. a large compar- 
ison sample with accurate spectroscopic redshifts and fully 
representative for the PACS fiux selection) at high z does 
not allow to properly quantify this statement. 
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Figure 7. Total IR Luminosity Function estimated by combining the data from the four PEP fields using the |Avni fc BahcaII| | |1980| l 
method (black filled circles). The grey filled area shows the uncertainty locus obtained by combining the Poissonian error with the 
photometric redshift uncertainty derived through Monte Carlo simulations. The black solid line represents our best fit to the PEP data 
in the diff'erent redshift bins, corresponding to the parameters report ed in Table [T] Other results from the literature are plotted for 
comparison (diagonal crosses connected by a grey dashed line, L LP of|Sanders et al.||2003[ magenta fille d stars, |Le Floc'h et al.||2005[ 
orange filled squares, [Chapman et al.|2005| orange dot-dash ed line,|Caputi et al.|2007| blue open trian gles, [Rodighiero et al.pOlOal green 
dashed line, [Magnelli et al. |2009| 2011; green open circles, [Magnelli et al.||2011[ l. Note that the Mag nelli et al.| l |2011| ) data correspond 
to slightly different redshift bins: 1.3<2<1.8 and 1.8<2<2.3, respectively. We have therefore plotted the data points corresponding to 
the first redshift bin in our 1.2<2<1.7 panel and the data point corresponding to the second redshift bin in both our 1.7<2:<2.0 and 
2.0<2<2.5 panels. The red open squares arc the total IR LPs derived by Vaccari et al. (in preparation) from the HerMES 250-/jm selected 
COSMOS sample. 



In Fig. [7] the total IR LFs obtained by combining the 
160-/xm selected samples with the Avni & Bahcall ( 1980 1 



technique is plotted and compared with other derivations 
available in the literature. The total IR LF of [Sanders eti 



al. (20031 is plotted as a local reference, in addition to the 
LFs of |Le Floc'h et al.] ([20051), |Rodighiero et aT] (|2010a|), 



Caputi et al. (20071, Magnelh et al. (20091 and Magnelh et 



al. (20111 in various redshift intervals. Globally, data from 



surveys at different wavelengths agree relatively well over 
the common z-range. No data for comparison are available 



at 0>2.5, apart from the IR LF of sub-mm galaxies from 



Chapman et al. ( 2005 1 and Wall et al. ( 2008 1 at 2~2.5, which 



represent reasonably well just the very bright end of the to- 
tal IR LF. Our derivation is the first at such high redshifts, 
especially in the 3<z^4.2 range. Note the good agreement 
between our PEP-based total IR LF and the HerMES-based 
one derived by Vaccari et al. (in preparation), shown by the 
red open squares in Fig. [7] Though consistent within the 
error-bars, in the highest redshift bin the HerMES LF is 
slightly higher than ours. Since the 250-^(m HerMES selec- 
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Figure 8. Total IR Luminosity Function estimated by combining the data from the four PEP fields using the |Avni fc Bahcall| | |l980| l 
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the itlcr (Poissonian) uncertainty regions at different redshifts. 



tion favours the detection of higher redshift sources than the 
PEP one, it is more hkely that the PEP LF in the higher-z 
bin is affected by some flux/redshift incompleteness rather 
than by the presence of low-z sources erroneously placed at 
high-2; by incorrect photometric redshift assignment (catas- 
trophic failures). The values of our total IR LF for each 
redshift and luminosity bin are reported in Table 6. 



3.4 Evolution 

In order to study the evolution of the total IR LF, we derive 
a parametric estimate of the luminosity function at different 
redshifts. For the shape of the LF we assume a modified- 
Schechter function (i.e. Saunders et al. 1990), where ^{L) is 
given by 
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behaving as a power law for L <^ L* and as a Gaussian 
in logZ/ for L ^ L* . The adopted LF parametric shape 
depends on 4 parameters {a, a, L* and $*), whose best 
fitting values and uncertainties have been found using a non- 
linear least squares fitting procedure. In detail, while in the 
first z-hm aU the parameters have been estimated, starting 
from the second z-bin, the values of a and cr have been 
frozen at the values found at lower redshift, leaving only L* 
and <I>* free to vary (see Table [t]). Note that, in the highest 
redshift bin (3.0<z<4.2) we are not able to constrain the LF 
break, while we are up to ~3. Therefore, the resuhs found 
at this redshift are affected by larger uncertainties than the 
z ^ 3 ones. However, although there is some degeneracy in 
the values of L* and "I>* at 3.0<z<4.2, the range of allowed 
value combinations still giving a reasonable fit to the three 



Figure 9. Evolution of L* and <&* as function of z (L*o<:(l -I- 
^■)3.55±0.10 2 < 1.85, <x(l + 2)i-62±0-5i at 2 > 1.85; <l>*(x(l -f 
2)-0-57±0.22 at 2<1.1, oc(l-H2)-3-92±0-34 at 2 > 1.1. 



observed data points (constraining the bright-end of the LF) 
is limited and do not significantly affect our results. 

In Fig. [8] we plot the total IR LFs at all redshifts with 
the ±l(j Poissonian uncertainty regions (different colours for 
different z-intervals). There is a clear luminosity evolution 
with redshift, at least up to z~3. The apparent "fall" of 
the 0>3 LF with respect to those at lower redshifts, if real, 
might imply a global negative evolution of the IR galax- 
ies and/or AGN starting at 2;>3. However, as mentioned 
above, we must point out that in the highest redshift bin 
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Table 7. Parameter values describing the curve fitted to the total IR LF 



redshift range 


a 


a 


logio(L*/LQ) 


logio(<I>*/Mpc"3 dcx^i) 


0.0<z<0.3 


1.15±0.07 


0.52±0.05 


10.12±0.16 


-2.29±0.06 


0.3<2<0.45 


1.2" 


0.5" 


10.41±0.03 


-2.31±0.03 


0.45<z<0.6 


1.2" 


0.5" 


10.55±0.03 


-2.35±0.05 


0.6<2<0.8 


1.2" 


0.5" 


10.71±0.03 


-2.35±0.06 


0.8<2<1.0 


1.2" 


0.5" 


10.97±0.04 


-2.40±0.05 


1.0<z<1.2 


1.2" 


0.5" 


11.13±0.04 


-2.43±0.04 


1.2<2<1.7 


1.2" 


0.5" 


11.37±0.03 


-2.70±0.04 


1.7<2<2.0 


1.2" 


0.5" 


11.50±0.03 


-3.00±0.03 


2.0<2<2.5 


1.2" 


0.5" 


11.60±0.03 


-3.01±0.11 


2.5<2<3.0 


1.2" 


0.5" 


11.92±0.08 


-3.27±0.18 


3.0<2<4.2 


1.2" 


0.5" 


11.90±0.16 


-3.74±0.30 



fixed value 



PACS data could be affected by incompleteness, since with 
increasing redshift (and for intrinsically "cold" sources) the 
true PACS fluxes might fall below the detection limit (i.e. 
faint sources should be missed even if completeness correc- 
tions are perfect), while, if luminous enough, they can still 
be detectable by SPIRE. This effect is expected to be more 
relevant in COSMOS, where SPIRE data are quite deep rel- 
atively to PACS, while it should not happen in GOODS-S, 
where PACS data are very deep compared to the SPIRE 
ones (which are limited by confusion). The high-« incom- 
pleteness of PACS surveys might also be emphasised by the 
redshift incompleteness of the COSMOS sample, that could 
affect the highest redshift bins more than the lower ones (al- 
though the redshift incompleteness seems to be independent 
Berta et al. 112011 1. However, we must point 



of redshift; see 



out that a decrease at 2^2.7-3 similar to that observed in 
our data, is also observed in the space density of X-ray (see 
Brusa et al.|2009 Civano et al.|20lT | and optically selected 
AGN ([Richards et al.||2006[ ), and of sub -mm galaxies JWalll 



|et al.|2008[ ), as weU as in the HerMES total IR LF - though 
less evident - (from 250-/im data; Vaccari et al., in prepara- 
tion: see Fig. [7|. Moreover, our result is in agreement with 
the recent finding of |Smit et al.| ( [2012[ ), that the character- 
istic value of the galaxy SFR exhibits a substantial, linear 
decrease as a function of redshift from 2;~2 to 2:~8. 
In Fig. |9] we show the values of L* and "I>* at different red- 
shifts, with the best curve (cx (l-l-z)*) fit to the data points. 
The values of the curve slopes and of the redshifts corre- 
sponding to evolutionary breaks have been chosen to be 
those which minimise the of the fit with two power-laws. 
We find a significant variation of L* with z, which increases 
as (1 + 2)3-55±o.io ^^1.85^ and as (1 + 2)i 62±o.5i 

tween 2~1.85 and 2~4. The variation of $* with redshift 
starts with a slow decrease as (1 + 2;)~''-5''='=0-22 2;~1.1, 
followed by a rapid decrease oc(l+2;)^^'^'^*'' '^'' at 2:> 1.1 and 
up to 2:~4. 

Previous estimations of the evolution of L* and "1>* (i.e. Ca- 



puti et al.|[2007| [Bethermin et al.|[20TT| and [Marsden et~ 



2011) discussed a decrease in the density of far-IR sources 
between 2=1 and 2=2. In particular, Bethermin et al. ( 2011 1 



and Marsden et al. (2011 1, by evolving a parameterised far- 



IR LF, explored the evolution required by the source counts 
in the parameter space. The results of these works (espe- 
cially those of [Bethermin et al.||20lT| showing a decrease 



of $* at 2>1 and a fiatter trend on the evolution of L* at 
2>2), are close to ours, though with the source counts only 
it was not possible to constrain the evolution of IR sources 
at 2>2. 

3.5 Evolution of the Different IR Populations 

The PEP survey, given its size and its coverage in lumi- 
nosity and redshift, allows us to go further in investigating 
the evolution of the total IR LF: it gives us the opportu- 
nity to study the evolution of the different galaxy classes 
that compose the global IR population. To investigate the 
different evolutionary paths of the various IR populations, 
we have computed the 1/Vmax LFs separately for the five 
galaxy classes defined by the SED-fitting analysis. In Fig. |10| 
we show the total IR LFs derived from the combined PEP 
samples for the different SED classes (coloured filled areas). 
The results of the fit to a modified Schechter function (see 
equation [3| for each population are overplotted on the data. 
In fact, by following the same procedure adopted for the 
global luminosity function, a parametric fit to the LFs at 
different redshifts has been performed also for the single 
populations. The a and a parameters, for each population, 
have been estimated at the redshift where the correspond- 
ing LF is best sampled (not necessarily at the lowest 2-bin 
as for the global LF). Subsequently, the values of a and a 
have been frozen at the values found in the "optimal" red- 
shift bin, leaving only L* and <I>* free to vary. In Fig. |11[ 
analogously to Fig. [oj we show the values of L* and <1>* at 
different redshifts for the different populations, with the best 
least square fitting curves (oc (I-I-2)'" ) overplotted. The best- 
fitting values of a, a, L* and in the first redshift bin, to- 
gether with the parameters describing the luminosity (fcL,i, 

fcL,2 and 2b, l: Oc(1-|-2)'''-'1 to 2=2b,L, 0c(1+2)'°l,2 at 2>2b,L) 

and density evolution {kp^\, fcp,2 and 2b, p: oc(l+2)'^'''i out to 
z=2b.p, oc(H-2)''f'^ at 2>2b,p) are reported in Table [s] 

A clear result of our analysis is that the evolution de- 
rived for the global IR LF is indeed a combination of dif- 
ferent evolutionary paths: the far-IR population does not 
evolve all together "as a whole" , as it is often assumed in the 
literature, but is composed by different galaxy classes evolv- 
ing differently and independently. As shown in Fig. |10| the 
normal spiral galaxy population dominates the luminosity 
function at low-2, from the local Universe up to 2~0.5. Mov- 
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Figure 10. Total IR Luminosity Function estimated with the 1/Vmax method by combining the data from the four PEP fields for the 
different populations (their ztlcr uncertainty regions are shown as coloured filled areas: green for spirals; cyan for starbursts; red for 
SF-AGN; magenta for AGN2; and blue for AGNl), compared to the total LF (ztlcr, grey filled area, same as in Fig.[7]|. The best-fit modified 
Schechter functions are also plotted, extrapolated to fainter and brighter luminosities than covered by the data (with the black curve 
being for the total LF and the same colours used for the filled areas as for the single populations). 



Table 8. Parameter values describing the curve fitted to the total IR LF of the different SED populations 



logio(L*/ 

L©) 
(0.0<z<0.3) 



logio(**/ 
Mpc~^ dex~l) 



2b, L 



V,2 



2b, p 



spiral l.OOitO.05 0.50±0.01 9.78±0.04 -2.12±0.01 4.49±0.15 

starburst 1.00±0.20 0.35±0.10 11.17±0.16 -4.46±0.06 1.96±0.13 

SF-AGM 1.20±0.02 0.40±0.10 10.80±0.02 -3.20±0.01 3.17±0.04 

AGN2 1.20±0.20 0.70±0.20 10.80±0.20 -5.14±0.17 1.41±0.33 

AGNl 1.40±0.30 0.70±0.20 10.50±0.20 -5.21±0.11 1.31±0.02 



0.00±0.46 



1.1 



-0.54±0.12 
3.79±0.21 
0.67±0.05 
2.65±0.32 
3.00±0.25 



-7.13±0.24 0.53 
-1.06±0.05 1.1 
-3.17±0.15 1.1 



ing to higher redshifts, the number density of galaxies with 
spiral galaxy SEDs sharply decreases, while their luminosity 
continues to increase, at least up to 2~1 (see the $* and L* 
parameter trends shown in Fig. We note that what we 
observe between z~0 and z^l for the spiral SED galax- 
ies is an increase of L* by a factor of ~5, and a decrease 
of $* by a factor of ~10. Since the two evolutions are not 
independent, the "total" evolutionary effect results from the 



combination of the two (as can be observed in the total IR 
luminosity density, see Section |4|. A way to derive the "to- 
tal" effect of evolution on a LF is to fix at a given volume 
density value and see how the luminosity corresponding to 
that value changes: indeed we find an increases by a factor 
of ~2.5 between z=Q and z=l for the spiral LF, in good 
agreement with previous results, either empirical (for mor- 
phologically classified disky galaxies; Scarlata et al.|2007 1 or 
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theoretical (from chemical evolution models of Milky Way- 
like galaxies; Colavitti, Matteucci fc Murante(2 008!). 
Over the whole redshift range 0.5—3, the "totall" luminos- 
ity function is dominated by the SF-AGN population. The 
number density of SF-AGN is nearly constant from the local 
Universe up to z~l-1.5, showing a slight decrease at higher 
redshifts, while their luminosities show positive evolution up 
to the highest redshifts (2:~3.5-4). From Fig. 10 we note a 
sort of bimodality in the SF-AGN LFs (at z < 0.45, where we 
are able to cover a larger luminosity range) . This bimodality 
is indeed to be ascribed to the crossing of two contributions: 
that of the SF-AGN(Spiral) population, responsible for the 
faint-end steepness of the LFs, and that of the SF-AGN (SB) 
population, dominating the bright-end of the SF-AGN LFs 
and declining at low Lm. (not reported in the figure). 
The starburst galaxy population never dominates. The red- 
shift range where we observe the highest contribution from 
the starburst galaxies is at z^l~2, while in the local Uni- 
verse their contribution is almost negligible (i.e. their 
parameter shows an opposite trend with respect to that of 
spiral galaxies see Fig. 



Ill 



The AGNl and AGN2 populations show a very similar evolu- 
tionary trend as a function of z, both in <1>* and L* . These 
powerful AGN populations dominate only the very bright 
end of the total IR LF, although their number densities and 
luminosities keep increasing from the local Universe up to 
the higher redshifts. At z>2.b the AGNl and AGN2 popula- 
tions become as important as the SF-AGN one, with the total 
IR LF of PACS-selected sources in the redshift range 2.5-4 
being totally dominated by objects containing an AGN. 



3.6 Total IR LF in Mass and Specific 
Star-Formation Rate bins 

3. 6. 1 Stellar masses and SFR from SED fitting 

The wealth of multi- wavelength data available in the cosmo- 
logical fields included in our work allow us to perform a de- 
tailed SED fitting of all sources, in order to derive their most 
relevant physical parameters (e.g. stellar masses). To derive 
stellar masses we have fitted the broad-band SEDs of our 



sources using a modified version of MAGPHYS ( Da Cunha et al 



2008 1 , which is a code describing the SEDs using a combina- 



tion of stellar light and emission from dust heated by stellar 
populations. In particular, the MAGPHYS software simultane- 
ously fits the broad-band UV-to-far IR observed SED of each 
object, ensuring an energy balance between the absorbed 
UV light and that re-emitted in the far-IR regime. The main 
assumptions are that the energy re-radiated by dust is equal 
to that absorbed, and that starlight is the only significant 



source of dust heating. We refer to Da Cunha et al. ( 2008 1 



for a thorough formal description of how galaxy SEDs are 
build. At each source's redshift, the code chooses among 
different combinations of star formation histories, metallic- 
ities and dust contents, associating a wide range of optical 
models to a wide range of infrared spectra and comparing 
to observed photometry, seeking for minimization. Each 
star formation (SF) history is parameterised in terms of an 
underlying continuous model with exponentially declining 
star formation rate (SFR), on top of which are superim- 
posed random bursts (see Da Cunha et al.|2008 Da Cunha 



|et al.||2010| ). We note that, although the MAGPHYS assump- 
tion of exponentially declining SFR might not be the best 
to reproduce the SFR history of 2>1.5 star- forming galax- 
ies (i.e. exponentially increasing or increasing SFR would be 
better choices, as widely discussed by |Maraston et aL]|2010| 
and |Reddy et aL||2012| ), in our specific case it does not af- 
fect the results. In fact, we do not use the MAGPHYS derived 
SFRs, but we compute them by integrating the best-fitting 
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z= 1.22 class = SF-AGN z= 0.38 class = AGN2 z= 1.25 class=AGN1 




Figure 12. Examples of PEP AGN SEDs fi tted by the MAGPHYS code ( |Da Cunha et al.|2008"l l modified by |Berta et al.| ( |2013| l to add an 
AGN component from tlie |Fritz et al.| | |2006[ l library. From left to right the result of the fit to a SF-AGN, an AGM2 and an AGNl are shown. 
The stellar component unattenuated by dust is shown in green, while the dust-attenuated spectrum with dust IR re-emission is shown 
in blue. The AGN contribution is shown in red, while the black curve is the total fitted spectrum, obtained from the sum of the blue 
and red components. Our data are represented by the black dots. 



SED (resulting from Le Phare). The models are distributed 
uniformly in metallicity between 0.2 and 2 times solar. Since 
the MAGPHYS code assumes that starlight is the only signif- 
icant source of dust heating, thus ignoring the presence of 



a possible AGN component, Berta et al. (20131 have devel- 



oped a modified version of the MAGPHYS code by adding a 
torus component to the modelled SED emission, combining 



the Da Cunha et al. ( 2008 1 original code with the Fritz et al. 
( 2006[ ) AGN torus library (see also |Feltre eraLl|2012[ ). The 



spectral fitting is performed by comparing the observed SED 
of our galaxies to every model in the generated library, at 
the corresponding redshift. A minimisation provides the 
quality of each fit. We must point out that the mass deriva- 
tion for unobscured AGN (i.e. AGNl) is a problematic issue, 
therefore the masses estimated for that class of objects are 
the most uncertain ones. One source of uncertainty in the 
mass measurement for AGNl is due to the fact that, in these 
objects, the UV part of the SED is likely dominated by the 
AGN rather than by the host galaxy. This may produce an 
underestimate of the mass, since, if the AGN contribution 
is not taken into account, the data can be fit by a bluer, 
younger and smaller mass object. On the other hand, the 
mid-IR part of the AGN SED is dominated by dust emis- 
sion from the dusty torus heated by the central black hole. 
If a proper decomposition into a stellar and a torus com- 
ponent is not performed, the use of a pure stellar template 
for estimating the mass from SED-fitting tends to reproduce 
the mid-IR emission with an older, redder and larger mass 
object (mass sometimes larger by a factor of 2 than those 
derived through a decomposition procedure; |Santini et al.| 
2012). These two effects lead to an increase in the uncer- 



tainty of the mass derivation, although they might somewhat 
compensate their effects for a large sample of objects. For 
this reason, we obtained measurements of the stellar masses 
of our objects containing an AGN by means of the specific 



decomposition technique developed by Berta et al. (20131, 



to separate stellar and nuclear emission components. Exam- 
ples of the results of this decomposition applied to SF-AGN, 
AGN2 and AGNl are shown in Fig. [l2] Masses of AGN esti- 
mated with the original MAGPHYS and with the |Berta et al.| 
( |2013[ ) code have been compared, showing very good agree- 
ment and small dispersion around the 1-1 relation. Similarly, 



as further check, we have also computed stellar masses with 
different code (Hyperz; BolzoncUa et al.|2000 1 and stellar li- 
brary (BC03, [Bru"zual fc Charlot|2003| instead of the CB07 
used as default by MAGPHYS), finding good agreement and no 
sistematics, too. 

We have derived the SFRs from the total IR lumi- 
nosities (estimated from the SED fitting described in Sec- 
tion 2.1 1 with the standard Kennicutt (19981 relation (con- 



verted to Chabrier IMF), after subtracting the AGN con- 
tribution to LiR. We note that the total IR luminosity in 
PEP sources is usually dominated by star formation, even 
in objects for which an AGN dominates the optical/near- 
IR/mid-IR part of the spectrum. 

3.6.2 LFs in different mass bins 

We compute the total IR LF for galaxies of different stel- 
lar masses: 8.5^1og(M/MQ)<10, 10^1og(M/MQ)<ll, and 
ll^log(M/M0)<12, by means of the standard 1/Vmax for- 
malism, and we show the results in Figs. [Ts] (total IR LF in 
different z-bins) and 14 (ratio between the LF in the mass 
intervals and the total IR LF). 

We have compared our results with the SFR function 



(SFR converted to IR luminosity using the Kennicutt ( 1998 1 



relation) of massive (log(M/M0>lO) galaxies derived by 
Fontanot et al.| ( |2012[ ) from the GOODS-MUSIC sample at 
redshift 0.4<2:<1.8. In the common redshift and luminosity 
range we find an excellent agreement with our total IR LF, 
which is dominated by sources with 10<log(M/MQ)<ll. At 
lower luminosities, not sampled by our data, the |Fontanot et| 
al. (20121 LFs are characterized by a double-peaked struc- 



ture, interpreted in terms of the well-known bimodality 
in the colour(SFR)-Mass diagram. As expected from the 
SFR-Mass relation, the knee (L*) of our IR LFs in differ- 
ent mass bins moves to higher luminosities with increasing 
masses (i.e. at 0.0<2:<0.3, log(L*/LQ) changes from ~9.5 
for log(M/M0)=8.5-lO sources, to ~11.3 for sources with 
log(M/MQ)>ll). 

The slope of the total IR LF in each mass bin is always 
similar to (or fiatter than) the "global" LF (total, includ- 
ing all the masses). The lower mass galaxies dominate at 
lower luminosities (log(L/LQ)<9), while the most massive 
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Figure 13. Contribution to the total IR LF from galaxies in three different mass intervals: 8.5<log(M/M0)<lO (green); 
lO<log(M/M 0)<ll (orange); and ll<log(M/M0)<12 (red). For comparison, the SFR function of log(M/MQ)>10 sources in the 
GOODS-S by |Fontanot et al.|2012| (converted to total IR LF) is showrn as black open squares. 



galaxies (log(A//M0)>ll) contribute only at higher lumi- 
nosities, even if they never dominate the LF. At all masses, 
the LF evolves with redshift, following the evolution of the 
"global" LF. Fig. |14| shows that the main contribution (>50 
per cent) to the total IR LF is due to intermediate-mass ob- 
jects (log(M/M0)=lO-ll) at all redshifts and luminosities, 
with their fraction remaining almost the same from z —0 
to z~4, simply shifting to higher luminosities. Lower mass 
objects (log(M/M0)=8.5-lO) contribute significantly only 
at log(I//L0)<lO, with their fraction just shifting to higher 
luminosities with redshift, but always being below 30 per 
cent at z>0.45 and log(L/L0)>ll. The contribution of the 
most massive objects (log(M/M0) = ll-12) increases with 
IR luminosity and redshift, becoming significant (>50 per 
cent) only at z>1.7 and log(L/L0)>12.5. Thus the bulk 
of the IR luminosity is produced by star-forming galaxies 
of mass around the characteristic mass M* of the Schecter 
mass function. 



3.7 Specific-SFR and the main sequence of 
star-forming galaxies 

Having computed stellar masses and SFRs for each source, 
we can check how the PACS selected sources populate the 
SFR-stellar mass plane and the so called main sequence 
(MS) of star-forming galaxies, as a function of redshift. This 
relation (i.e. SFR versus stellar mass) has been shown to be 
quite tight in the local Universe (Peng et al. 2010 2011) 
and well established at redshift 



A (Elbaz et al 

and up t o z~2 ([Daddi et al.|[2007l [Rodighiero et al 



20071 



20111 



and ( Magdis et al.|20"l0 l, with normalisation scaling as 



~(1+^)^- 


" out to «~2, as shown by 


Sargent et al. 


(2012| (see 


also|Elbaz et al.||2007| 


Rodighiero et al.||2010b| | 


-"annella et| 


al. 2009 


Karim et al. 


2011 


I. At 2~2 and ~1.5 we assume 



a slope of 0.79 for the MS in the SFR versus stellar-mass 
plane (according to [Rodighiero et al.]|2011| and [Sargent "et| 
al.||2012|), while at z~l we assume a slope of 0.9, as found 



by (Elbaz et"al] ( [2007| l, and we limit our investigation to the 
redshift range 0.8<z<2.2. 

By combining UV and far-IR data, [Rodighiero et aL] 
( [2011[ ) re-evaluated the locus of the MS at z~2, showing that 
objects lying a factor of 4 above the MS (in SFR) can be con- 
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Figure 14. Ratio between the IR LF of galaxies in three different 
mass intervals: 8.5<log(M/Mo)<10 (top); lO<log(M/M0)<ll 
{middle); and llsglog(M/MQ)<12 {bottom), as shown in Fig. [13] 
and the total IR LF, plotted in all the different redshift intervals 
considered in this study, from 2^0 to Z'~4. The different colours 
represent the different redshifts, as explained in the plot. 



sidered as outliers with respect to the average locus where 
smoothly star-forming galaxies spend most of their lives in a 
secular and steady regime. In that work, off-sequence sources 
(characterized by very high specific-SFRs) are assumed to 
be in a starburst mode, and are found to contribute only 
2 per cent of mass-selected star-forming galaxies and to ac- 
count for only 10 per cent of the cosmic SFR density at 
i Rodighiero et al.|20H" |. 

In order to check what kind of objects we could classify 
as on- and off-MS sources in our IR sample compared to 
previous findings, based either on IR or on optical surveys 



(e.g. Rodighiero et al. 2011 Sargent et al. 20121, we have 
splitted our sample into off-MS and on-MS. For consistency 
with previous studies we have applied the same criterion as 



Rodighiero et al. (20111 (0.6 dex above the MS) over the 



whole 0.8<z<2.2 redshift range, by using as a reference MS 
the relation found by Rodighiero et al. (2011 1 at z^2, scaled 
as described above at z 
let al.l ( pOOT] ) at 



1.5, and the relation found by |Elbaz"l 
In Fig. |15[ we show the SFR versus 



stellar mass distributions in three redshift bins (0.8<2;<1.25, 
1.25<2:<1.8 and 1.8<2<2.2), for the PACS sources included 
in the computation of the luminosity functions presented in 
this work. The colour code marks the different SED-classes 
to which each source belongs. We also report the typical loci 
of the MS at the various redshifts (scaling as (l-|-z)^'*, as 
mentioned above). Details are given in the caption of the 
Figure. The typical far-IR selection bias (PACS-Herschel in 
this case) appears as an approximate horizontal SFR cut 
2011[ |Wuyts et aI1|2011[ ), shown as thin 



(Rodighiero et al. 



dotted line in Fig. 15 We note that the trends of mid/far-IR 
SEDs with offset from the main sequence observed in Fig |15| 
(and widely discussed in Section [5| are in good agreement 
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Figure 15. SFR versus stellar mass for our PEP 160-^m 
sources (green, spiral; cyan, starburst; red, SF-AGN; magenta, 
AGN2; blue, AGNl), in three redshift bins (from left to right): 
0.8<2;<1.25; 1.25<2<1.8; and 1.8<2<2.2. The relation known as 
main sequence is plotted as a solid line (from |Elbaz et al.||2007| 
in the lower redshift bin, rescaled as {1+z)^'^ in the central bin 
and from [Rodighiero et al.||2011| at 2~2), while the dashed-line 
shows the same relation 0.6 dex higher, indicating the separation 
between MS and above MS objects adopted by [Rodighiero et al.| 
| |2011[ l. The horizontal dotted lines show the nominal SFR lim- 
its of the GOODS-S (lower) and COSMOS (upper) samples in 
the different redshift intervals. The orange open squares and the 
dark-green open diamonds mark the two sub-classes of SF-AGN 
galaxies: the SF-AGN (SB) and SF-AGN (Spiral) respectively. 



with the results of Elbaz et al 
( [2012 1 



(20111 and Nordon et al. 



With the selection based on [Rodighiero et aX\ \2Q11\ 



and overplotted in Fig. 15 (sources qualify as "off-MS" if 



they lie more than 0.6 dex above the observed SFR-stellar 
mass relation) , we can compute the contribution of off- 
MS (also called "starburst" in the literature) and on-MS 
("steady star-formers" in the literature) sources to the to- 
tal IR LFs. This is presented in Fig. |16| where the total IR 
LFs of on- and off-MS sources have been computed indepen- 
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Figure 16. Contribution of MS (ztltr uncertainty region: pink shaded area) and ofF-MS (itlfT uncertainty region: yellow shaded area) 
galaxies to the total IR LF (black filled dots and area) in three difi'erent redshift bins. For comparison, the recent esimates of |Sargent et| 
|al.| | |2012] | in similar ^-bins are shown as grey (total LF, in background), light grey (MS LF) and dark grey (off'-MS LF) filled regions. 



dently for the three redshift bins. The pink and yellow filled 
areas correspond to the ila uncertainty regions of the to- 
tal IR LFs estimated for the on-MS and off-MS populations, 
respectively, while the black filled area marks the global pop- 
ulation. Our results are compared with the recent estimates 



by Sargent et al. (20121, a non-parametric approach that is 
based on three basic observables: the redshift evolution of 
the stellar mass function for star-forming galaxies; the evo- 
lution of the sSFR of MS galaxies; and a double-Gaussian 
decomposition of the sSFR distribution at fixed stellar mass 
into a contribution (assumed redshift- and mass-invariant) 
from MS and off-MS (i.e. starburst) activity. The evolution 
of the two populations found both in our data and the |Sar-| 
gent et al. (|2012l) model are very similar. Both data and 
estimates indicate that the bright-end of the total IR LF is 
dominated by off-MS sources. However, although consistent 
within the uncertainties, the relative contribution of off-MS 



sources seems to be stronger in the Sargent et al. (20121 



model than observed in the present computation, where we 
find that the bright-ends of the PEP IR LFs are more sim- 
ilarly populated by MS and off-MS sources (especially at 
2;~2). This difference can be at least partly ascribed to the 
sharp cut we apply to separate MS from off-MS sources, 



while Sargent et al. (20121 model the off- and on-MS popu- 



lations with two continuous log-normal distributions centred 
at 0.6 dex above the MS and exactly on the MS respectively, 
of which the one describing the "starburst" population has 
wings that extend into our on-MS selection region (hence 
attributing more sources to the starburst category than are 
selected in our off-MS class). A better agreement between 



our data and the estimates of Sargent et al. ( 2012 I is found at 



the faint-end of the LFs that appear to be completely domi- 
nated by the "normal" MS galaxies at all redshifts (although 
the total and relative contributions at log(Lia/L0)<12 are 
slightly lower in the data than predicted by [Sargent et al.| 
2012 1. Good agreement is also found with respect to the evo- 



lution of the cross-over luminosity (i.e. where the contribu- 
tions from on- and off-MS sources are equal); in the model of 
Sargent et al. (20121 the cross-over luminosity simply shifts 



luminosity and density evolution considered). Note that the 
model assumptions rely on results from different surveys, se- 
lected at different wavelengths, complete in mass and with 
good sampling of the MS. On the other hand, our selection 
is in SFR and our sources do not follow any clear sequence 
in stellar mass-SFR plane, because, except at the highest 
masses, the data are not deep enough to reach well into the 
main sequence. These different selection effects are likely to 
lead to some differences between our LFs and the estimates 



of Sargent et al. (20121 



To quantify the relative contribution of the two popu- 
lations, our observed data have been fitted with a modified 
Schechter function, in order to integrate them and compute 
their comoving number and luminosity densities as functions 
of redshift (see next Section). 



4 NUMBER DENSITY AND IR LUMINOSITY 
DENSITY 

We derive the evolution of the comoving number and lu- 



minosity density (total, see Fig 18 1 of the PEP sources, ei- 



ther belonging to the different SED classes (Fig 18 left), 



to higher luminosities and lower densities (according to the 



to the on- and off-MS categories (middle) and to the dif- 
ferent mass intervals (right), by integrating the total IR 
LF in the different redshift bins from z~0 to z^4. To com- 
pute the number (and IR luminosity) density, we integrate 
the Schechter functions that best reproduce the different 
populations/mass/sSFR-classes, down to log(I//L0)=8. We 
note that here we consider lower limits the number and lu- 
minosity densities at 3.0<z<4.2, since our LF estimate in 
that redshift bin is likely to be incomplete, as discussed in 
Section fS.SI We find that the number density of the whole IR 
population is nearly constant in the 2=0-1.2 redshift range 
(slightly increasing from z~0 to 2~0.5), decreasing at z>1.2 
(see top panels of Fig. |18[ ). When decomposing the number 
density according to the different SED classes, we observe 
that normal spiral galaxies dominate the local density, with 
a smaller contribution also from the SF-AGN population and 
a negligible one from starburst, AGNl and AGN2. The space 
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density of spiral galaxies decreases rapidly at while 
that of SF-AGN stays nearly constant at 0.5 S z S 2.5, largely 
dominating in that redshift range. Starburst galaxies never 
dominate, while the number density of the bright AGN (both 
AGNl and AGN2) increases with redshift, from ~10~'' Mpc~^ 
at 2:~0 to ~l-2xl0"^ Mpc"^ at z~3. At higher redshifts 
the AGN population largely dominates the number density. 
If the overall contribution to the IR luminosity density 
(Pir) from the AGN components of galaxies is small, pm. 
can be considered as a proxy of the SFR density (psfr). 
As a further check, we have therefore studied the evolu- 
tion of the SF-AGN population (which dominates the distri- 
bution of sources) by dividing this class into SF-AGN (SB) 
and SF-AGN (Spiral) sub-classes and studying their evolu- 
tion separately. Indeed, we have found different evolution- 
ary paths for the two populations, the former dominat- 
ing at higher redshifts and showing a behaviour similar to 
that of AGN-dominated sources (e.g. AGNl and AGN2), the 
latter dominating at intermediate redshifts (between 
and 2), rising sharply from 2~2 toward the lower redshifts 
and decreasing, while the spiral population rises at z^l. 
These evolutionary trends, in terms of number and lumi- 
nosity density, have been reported in Fig. [18] as orange 
dot-dot-dot-dashed (SF-AGN (SB)) and dark-green dashed 
(SF-AGN(Spiral)) curves. 

Galaxies following the SFR-mass relation are always domi- 
nant over the off-MS population, at all redshifts (although 
their space density decreases with increasing z, as well as 
the "global" number density), while the number density of 
the latter population remains nearly constant between 2;~0.8 
and z~2.2. 

In all the mass bins, the trends with redshift of the galaxy 
number densities are similar to the "global" one, decreasing 
at higher redshifts, although with slightly different slopes 
for the different mass intervals. The number densities of 
low mass galaxies (8.5<log(Af/MQ)<10), reported in the 
top right panel of Fig. |18[ have been computed by inte- 
grating the best-fitting modified Schechter function only 
to z~2, since data were not enough to derive reliable fits 
at higher redshifts. To this redshift, these sources outnum- 
ber the higher mass ones, although they fall steeply above 
2~1, when they reach about the same volume density of 
higher mass gala:xies (lO<log(M/M0)<ll). Massive objects 
(log(M/MQ)>ll) never dominate (always below 5 per cent) 
the total number density. 

The total IR LF allows a direct estimate of the total 
comoving IR luminosity density (pm) as a function of z, 
which is a crucial tool for understanding galaxy formation 
and evolution. Although piR, can be converted to a SFR 
density (psfr) under the assumption that the SFR and Lir 



quantities are connected by the Kennicutt ( 1998 1 relation, 
before doing that we must be sure that the total IR lu- 
minosity is produced uniquely by star-formation, without 
contamination from an AGN. The SED decomposition and 
separation into AGN and SF contributions show a negligible 
contribution to Lir (<10 per cent) from the AGN in most 
of the SF-AGN, and a SF component dominating the far- 
IR even in the majority of more powerful AGN (AGNl and 
AGN2). Here we prefer to speak in terms of piR rather than 
of psFR, since, especially at high redshift - where the AGN- 
dominated sources are more numerous - the conversion of 
Pir could represent only an upper limit to psfr. Note, how- 



o 

CL 




Figure 17. Redshift evolution of the total IR luminosity den- 
sity (piR, obtained by integrating the Schechter functions that 
best reproduce the total IR LF down to log(L/L0)=8) to z=4. 
The results of integrating the best-fitting curve for our observed 
total IR LF in each z-bin are shown as black filled circles (the 
grey filled area is the itia- uncertainty locus) and compared with 
estimates from previous mid-IR surveys (magenta filled area, |Le| 
[Floc'h et al.|2005 | orange filled triangles, fCaputi et al.|2007| blue 
open triangles, j Rodighiero et al. 2010a| and green open circles, 
[Magnelli et al.|2011 ). The upward pointing arrow in the highest-z 
bin means that, due to the large fraction of photometric redshifts 
and the fact that the PEP selection might miss high-z sources, 
our 3.0<2<4.2 pjji estimate is likely to be a lower limit. 



ever, that since this population is never dominant in our IR 
survey, we do not expect that contamination related to ac- 
cretion activity occurring in these objects (mainly at high-z) 
can significantly affect the results in terms of psfr. 

In Fig. [it] we show piR estimated from our total IR LF 
and compare it with results obtained from previous IR sur- 



veys ( 


Le Floc'h et al.||2005 


Caputi et al.||2007 Rodighiero 


et al. 


2010a MagnelH et al. 


20111. In the common redshift 



intervals (O<0< 2-2.5), we find very close agreement with 
previous results based on IR data, especially with the |Mag^ 



nelli et al. ( 2011 1 derivation. As well as previous findings, piR 
from PEP shows the rapid rise from 2~0 to 2~1, followed 
by a flattening at higher redshifts. The indications from our 
survey are that the intermediate redshift flattening is fol- 
lowed by a high redshift decline, which starts around z~3. 
From our data, piR evolves as (l-t-z)^ *^*"'^ up to z^l.l, as 
{l+z)-°-'^^°-^ from 2-1.1 to 2-2.8, then as (l+2)-8-°±°-^ 
up to 2—4. 

In the bottom panels of Fig.[T8]we plot the different con- 
tributions to piR from the different SED populations (left), 
from the on- and off-MS sources (middle) and from the dif- 
ferent mass intervals. We notice a predominance of spiral- 
SED galaxies only at low redshifts (2<0.5-0.6), when SF-AGN 
begin to dominate piR up to 2—2.5. The starburst SED 
galaxies are never the prevalent population, although their 
contribution to piR increases rapidly from the local Universe 
to 2—1, then keeps nearly constant to 2—2.5, to decrease 
at higher redshifts. The SF-AGN(SB) and SF-AGN(Spiral) 
contributions to pm show opposite trends, with the former 
sharply increasing towards the higher redshifts (dominat- 
ing at 2>2), and the latter prevailing between 2—1 and —2, 
then dropping at higher redshifts. AGNl and AGN2 start dom- 
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Figure 18. Top: Evolution of the comoving number density of PEP sources up to (black filled circles with error-bars within the 
itlcr uncertainty region, represented by the grey filled area). The upward pointing arrow in the highest-z bin means that our 3.0<2<4.2 
estimates are likely to be lower limits. Bottom: Redshift evolution of the total IR luminosity density to z=4. To compute the number 
and IR luminosity density, we integrate the Schechter functions that best reproduce the different populations/mass/sSFR-classes, down 
to log(L/LQ)=8. The black filled circles and the grey dashed area in all the three panels represent our PEP derived pjp^ and its itlo" 
uncertainty region, as shown in Fig. |17[ In the left panels we show the number (top) and luminosity density (bottom) of the different 
IR populations (green filled area, spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; and blue, AGNl). The contribution of SF-AGN 
sources, sub-divided on the basis of their SED resemblance to spiral or starburst templates, are shown by the dark-green dashed 
(SF-AGN (Spiral)) and orange dot-dot-dot-dashed (SF-AGN(SB)) lines, respectively. In the middle panels we show the uncertainty regions 
of the relative contribution in number and luminosity density of the sources on- and off- the SFR-stellar mass MS jElbaz et al.[|2007[ 
[Daddi et al.|2007| , as pink and yellow filled areas, respectively. Our derivations have been compared to those of |Bethermin et al.| | |2012| 
for on- and off-MS sources (converted from psFR to PIR using the |Kennicutt|Pl998| relation), which are represented by the purple and 
orange dashed lines respectively. In the right panels we show the relative contribution to the number and luminosity density of sources 
with different masses (green, 8.5<log(M/MQ)<10; orange, 10<log(M/MQ)<ll; and red, ll<log(M/MQ)<12). For comparison, in the 
bottom right panel we plot also the results of |Santini et al.| ( |2009^ in the GOODS-S field in similar mass intervals (light-orange triangles 
and dashed line: 9.77<log(M/M0)<lO.77; pink triangles and dashed line: log(M/M0)>lO.77). 



inating the IR luminosity density at z it 2.5, with their pm. 
always rising from z^O, then remaining almost constant (or 
slightly decreasing) towards the higher redshifts. Note how 
these two populations of AGN evolve similarly, as an indi- 
cation of their same intrinsic nature and of the absence of 
any significant bias in the far-IR selection. 

The contributions to piR, of MS and ofT-MS sources 
(pink and yellow filled regions, respectively, in the 
bottom middle panel of Fig. |18[ ) stay nearly constant be- 
tween 2~0.8 and 2~2.2. In particular, the ofF-MS sources 
contribution stays around 20 per cent of the total piR over 
the whole 0.8<2<2.2 redshift interval, showing no signifi- 
cant signs of increase (or decrease). Our results strengthen 



the role of MS sources (pink shaded regions) in the build 
up of the stellar mass in galaxies, at all cosmic epochs, with 
evidence that their role is even increasing from z^2 to z~l: 
their number density changes by a factor of ~2, while their 
luminosity/ SFR density remains nearly constant. The im- 
portance of the off-MS sources does not show any significant 
signs of decreasing at lower z, their relative (with respect to 
the total) number (luminosity) density passing from ~9 per 
cent (22 per cent) at z~2 to ~6 per cent (19 per cent) at 
These fractions are relatively different from those found 



by Rodighiero et al. (2011 1 (off-MS galaxies represent only 2 



per cent of mass-selected star-forming galaxies and account 
for only 10 per cent of the cosmic SFR density at 2:~2). 
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However, we must note that Rodighiero et al. (20111 for 



their analysis combined far-IR-selected (i.e., SFR-selected) 
and near-IR-selected (i.e., M*-selected) star-forming sam- 
ples, well defining the main sequence, while with our data 
(SFR-selected only), we are not able to observe any correla- 
tion between SFR and stellar mass (see Fig. 15 I and barely 
detect MS objects at z~2. 

Our results have been compared to the SFR densities (con- 
verted to piR using the Kennicutt|1998 relation) for on- and 
off-MS sources based on the Sargent et al. (20121 model re- 



cently derived by Bethermin et al. (20121. These are shown 



in Fig. [18] as purple and orange dashed lines respectively. 
The predicted off-MS psfr agrees well with our estimate in 
the common redshift range, while the predicted one for MS 
sources is higher than that derived from our data (especially 
at 2;~1.5). As already discussed regarding the comparison 
with the Sargent et al. (20121 model, this discrepancy is 



likely to be ascribed to our difficulty to extrapolate the MS 
to low SFR values and to the different selection criteria used 
to separate MS from off-MS sources. 

In the bottom right panels of Fig. |18| we show the con- 
tribution of the different mass populations to the luminos- 
ity density as a function of redshift. Although we detect 
a similar steep increase of piR versus redshift at zUl in 
both low and high mass galaxies, the evolution in piR of 
galaxies with different masses is very different, reflecting 
the downsizing scenario, with piR peaking at higher red- 
shift with increasing mass. Indeed, the IR luminosity density 
of intermediate-mass objects (log(M/MQ) = 10-ll) always 
dominates, increasing up to z~l, then remaining nearly con- 
stant at higher redshifts (at least up to z~2.8). The IR lumi- 
nosity density of most massive objects increases even more 
rapidly with redshift (at z=2 it was ~5 times higher than 
today) and continues to grow up to z=3, where their contri- 
bution to piR is ~30 per cent of the total and close to that of 
intermediate mass objects (which contribute ~60 per cent 



at z>3). We compare our results with those of Santini et 



al. (20091, plotted in the figure as thick dashed lines (light 
orange: 9.77<log(M/MQ)<10.77; pink: log(M/MQ)>10.77; 
with a Chabrier IMF used to determine our masses), show- 
ing very similar trends and values for both intermediate- and 
higher-mass galaxies. Our analysis of high mass galaxies ex- 
tends up to z~4, finding that for the most massive galaxies 
piR continues to rise even at z^2, with an apparent peak 
at 2 = 3. This result confirms that the formation epoch of 
galaxies proceeded from high- to low-mass systems. 

The values of piR in the different redshift intervals, ei- 
ther the total ones or the contributions from the different 
classes (SED, mass, sSFR), are reported in Table 9. 



5 DISCUSSION 

In the previous sections we have discussed the different evo- 
lutionary behaviour of different classes of sources, either di- 
vided by SED-type, mass or sSFR. In this section we will try 
to understand "who is who" , discussing which populations 
are mainly on- or off-MS, which have the larger (smaller) 
masses, and how and if the relative contributions of these 
populations vary with redshift. 

From Fig. [15] we note that the off-MS sources are dom- 
inated by galaxies with AGN-typc SEDs. In the lower red- 



shift bin considered, MS sources are mainly spirals and 
SF-AGN(Spiral) , with a small tail of these populations con- 
tributing also to the off-MS. AGNl and AGN2 are prevalently 
0.6 dex above the MS (off-MS). Although the mass estimate 
of AGNl suffers of the largest uncertainties, our result sug- 
gests that the off-MS population is largely constituted by 
sources containing an AGN, with a higher fraction of AGN- 
dominated objects at higher z. This should indicate that the 
major merger episode likely associated with what is consid- 
ered a different mode of star-formation (e.g. Wuyts et al. 
|2011| finds that off-MS galaxies are very compact from their 
Sersic index, and quite likely mergers), could also trigger 
an intense AGN activity, whose presence strongly influences 
the SFR within the host galaxy. While AGN show preferen- 
tially high sSFRs, in our sample they do not seem to prefer 
higher- mass systems (see Fig. 191. However, this may again 
be due to our SFR selection, so we miss normal main se- 



quence galaxies especially at low mass (cf. Rodighiero et al 



20111 Therefore, our AGN are still in an intense SF phase, 
where they are still increasing their stellar mass by actively 
forming stars. Due to the far-IR selection, in fact, we miss 
the further phase of quiescent, more massive population, 
with the SF totally quenched and all the stellar mass al- 
ready formed and in place. Moreover, our result is not in 
conflict with the consensus in the literature (i.e. that X-ray 
AGN above a certain X-ray luminosity do prefer massive 
galaxies), since our classification is not in AGN luminosity 
cut, but more a cut in Lagn/ isF ratio. The recent result 
of Mullaney et al. (20121, that Lx/iiR (which is a proxy of 
^AGN / Lqy) is on average almost independent of stellar mass, 
implies that our selected AGN are not necessarily hosted by 
very massive galaxies, as indeed we find. 



Our results seem to confirm those of Santini et al. (20121, 
finding evidence of a higher average SF activity in AGN 
hosts with respect to inactive galaxies and a more pro- 
nounced level of SF enhancement in the hosts of luminous 
AGN (i.e. our AGNl and AGN2). Most of the starburst galax- 
ies (i.e. with SEDs fitted by local "starburst" templates) 
are classified on-MS at any redshifts, though they are al- 
ways in the region between the MS and the threshold. The 
starburst population seems to enhance its sSFR from z~2 
to z^l, occupying the region around the MS at higher z and 
shifting to higher SFRs (peaking at ~2x the MS) at lower 
z. The bulk of the spiral population remains around the 
MS at all redshifts, although they almost disappear from 
our sample at z~2. From Fig. [l5] we also note that, while 
the SF-AGN population as a whole occupies both the loci 
of MS and off-MS sources, when divided into SF-AGN(SB) 
and SF-AGN (Spiral) , it shows a very clear segregation (sug- 
gesting a different mode of star-formation for the two sub- 
classes). In fact, SF-AGN(SB) are mainly concentrated in the 
off-MS region, while the bulk of the SF-AGN(Spiral) is on- 
MS. 



If we mix together all the "ingredients" presented above, 
we can try to give a global interpretation to our results in 
terms of galaxy evolution. In agreement with other Herschel 



findings (i.e. Shao et al. 


2010] |Lapi et al. 


2011 


[Santini et 


al.||2012| IRosario et al.| 


2012| |Page et al.||2012| |Mullaney 


et al.]|2012l, we propose 


the following twofold evolutionary 



scenario for galaxies and AGN (sketched in a cartoon in 



Fig. 201: 
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Figure 19. Mass distribution of the PEP 160-fim sources in all the four PEP fields, with different colours corresponding to different 
SED-classes (green, spiral; cyan, starburst; red, SF-AGN; magenta, AGN2; and blue, AGNl). The top panel shows the mass distribution 
of the different fR populations at all redshifts, while the bottom panels report the mass distribution in three representative redshift bins 
(0<2<0.3, 0.8<2<1.0 and 2.0<z<2.5). 



• On the one hand, we observe that the sources with 
AGN-dominated SEDs (either AGNl or AGN2) and those with 
a starburst-hke SED, but containing a non-dominant AGN 
(SF-AGN (SB)) have a peak in number and luminosity den- 
sity at z~2-2.5, dominating at higher redshifts and rapidly 
decreasing at lower z. The evolution of AGNl and AGN2, 
both in luminosity and in density, is very similar, suggesting 
the same nature for type 1 and 2 AGN, and the unbiased 
power of observations in the far-IR band towards orienta- 
tion/obscuration (affecting both optical and X-ray observa- 
tions). While dominating the mid-IR part of the SED of 
these objects, the AGN is not able to explain the high ob- 
served far-IR emission, which is mostly powered by star- 
formation (as for the SF-AGN (SB) population, where the 
AGN never dominates the energetic output, probably due 
to dust-obscuration). The hosts of these AGN appear to 
form stars in a very efficient way, placing a large fraction of 



them above the known stellar mass-SFR MS (see Fig. 15 1. 
AGNl, AGN2 and SF-AGN (SB) are likely the progenitors of the 
elliptical galaxies observed nowadays in optical and near-IR 



surveys, forming through an intense burst of star-formation 
occurring during major mergers or in dense nuclear star- 



forming regions (Granato et al. 2001 Daddi et al. 2010 



Wuyts et al. 20111, then followed by a phase of nuclear 



activity during which their SMBHs grow (i.e. Hopkins et 
2008b: 



al 



2008a 



Lapi et al. 20111. It is generally agreed 



that the SMBHs and their host galaxies are tightly related, 
with major-mergers being considered the likely process re- 
sponsible for transporting large amount of gas towards the 
centre of the merging system, feeding the SMBH and trig- 
gering the SF activity. After the intense starburst phase. 



the AGN are believed to suppress the SF (i.e. Di Matteo 
et al. 20051, so that the remnant quickly evolves to a red 



massive spheroid. This picture is strongly supported by the 
recent Herschel results from PEP ( Rosario et al.|20l2 \ , Her- 



MES ([Page et al.|2012[ [Harrison et al.|2012[ ) and H-ATLAS 
( |Lapi et al.||201i r Surveys. The latter work, in particular, 
has shown that the bright-end of the IR LFs and counts 
at high redshift (>1.5) are consistent with the picture (e.g. 
[Granato et al.|[2001[ ) predicting the presence of a popula- 
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Figure 20. A cartoon showing a possible evolutionary scenario involving AGN and star-forming galaxies and leading to the formation, 
on the one hand, of local elliptical galaxies, on the other hand, of local spiral galaxies. Top: a strong starburst with a growing BH inside 
(SF-AGN(SB)) transforms into an AGN (either AGNl or AGN2, depending on orientation) when the BH mass and AGN luminosity have 
grown enough and the star formation is likely to be quenched by feedback processes. Thereafter the galaxy evolves passively toward a 
local elliptical. Bottom: the initial moderate starburst (starburst, lasting a typical time of ~few 10* yrs) transforms in a less intense 
starburst also fueling a low-luminosity AGN (SF-AGN (Spiral)), which starts to reveal itself when the starburst activity is fading. Once 
the AGN is triggered, it heats the dust resulting in an increase in the mid-infrared luminosity (flattening of the SED between 3 and 
8^m). The SF-AGN (Spiral) system could last about 7x longer than the pure starburst phase before becoming a steady spiral galaxy 
(spiral). 



tion of strongly-obscured, star-forming galaxies with SED 
appreciably different from those of the local starbursts. X- 



ray and mid-IR spectroscopic observations (e.g. Alexander 
erar||2005| [Vliante et al][2007l ) of sub -mm selected sources 



have revealed in many of them the presence of a growing 
SMBH, powering an obscured AGNs. In this framework, the 
most likely evolutionary path envisages first a SF-AGN (SB) 
phase (star- forming galaxy with a growing SMBH inside), 
then an AGN-dominated phase, and finally the formation of 
an elliptical galaxy in passive evolution. 

• On the other hand, we observe that a significant frac- 
tion of our IR selected sources is constituted by moderately 
star-forming galaxies characterised by an SED similar to 
that of spiral galaxies, but also suggesting the presence of a 
low-luminosity AGN (SF-AGN(Spiral) ), best-fitted by local 
Seyfert 1.8/2 templates. The bulk of these objects and their 
principal contribution to piR are at intermediate redshifts 
(lS:zJ:2), while they decrease between 2~1 and z=0, as 
the spiral population rises. Most of the SF-AGN (Spiral)-, 
spiral- and starburst-SED galaxies occupy the region 



around the SFR-stellar mass MS at any redshifts, suggesting 
a steady mode of star-formation rather than a "starburst" 
one for these three populations, whose evolution is likely 
connected. Indeed, MuUaney et al. (20121 have recently 



shown that at least up to z~2, SMBHs have grown together 
with their host galaxies in star-forming galaxies, irrespec- 
tive of host galaxy mass and triggering mechanism. Given 



the number densities (Fig. 18 1, evolutionary trends (Fig. 11 1, 
SFRs and masses (Figs. |15|and[T9l) of the SF-AGN(Spiral) , 
starburst and spiral populations, we suggest that these 
three classes of objects might constitute different phases in 
the life of a galaxy undergoing secular evolution. The gas 
in moderate starburst galaxies, undergoing a burst of en- 
hanced SF due either to gravitational interactions or disk in- 
stabilities (typical burst duration time of the order of a few 
10^ yr), might also fuel a low-luminosity AGN, which starts 
to reveal itself when the starburst activity is fading. Given 
the relatively high stellar masses found for the bulk of the 
SF-AGN(Spiral) (log(M/M0 
almost constant A/bh/M 



10-11, see Fig. 15 1, and the 
ratio (of 



-1-2x10"'') recently 
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suggested by Mullaney et al. (20121 for all the 0.5<2;<2.5 



star-forming galaxies, they are likely to contain relatively 
massive BHs (Mbh~10^-10* M©). They can therefore have 
either low values of their radiative efficiency or low values of 
their accretion mass rate m (m=M/MEdd<0.01). AGN with 
low radiative efficiency or low accretion mass rates are gen- 
erally called radiatively inefficient accretion flows (RlAFs - 
|Nara yan"1995-, which include also the advection-dominated 
accretion flow, ADAF - jNarayan 19941. Low m AGN are 
more difficult to detect, since they often are less luminous 
than their hosts. Because of this, the nuclear emission is di- 
luted by the host galaxy's emission and many of them are 
likely to be classified as "normal galaxies" in most surveys if 



the AGN luminosity is less than that of the host (e.g. Hop- 
kins et al.|2009 1. Given the relative number densities of the 



starburst and SF-AGN (Spiral) populations at l<z<2, we 
hypothesize a typical duration of the SF-AGN (Spiral) phase 
about 7 times longer than the starburst one (typical burst 
duration ~10* yr). Then, after a typical time of ~7xl0* 
yr, the AGN activity stops and these objects, whose num- 
ber density decreases at 0< 1, are likely to become steady 
spiral galaxies (rapidly increasing between 2~1 and z=0) 
at lower redshifts. 



6 CONCLUSIONS 

We have used the 70-, 100-, 160-, 250-, 350- and 500-^im 
data from the cosmological guaranteed time Herschel sur- 
veys, PEP and HerMES, in the GOODS-S and -N, ECDFS 
and COSMOS, to characterise the evolution of the IR lu- 
minosity function and luminosity density of PACS selected 
sources across the redshift range Q^z^A. Evolution is well 
constrained by our data up to 2~3, strong hints of evolution 
are derived at 3<2<4. In the present work we have: 

(i) completely characterised the multi-wavelength SEDs 
of the PEP sources by performing a detailed SED-fitting 
analysis and comparison with known template library of 
IR populations. Sources have been classified, based on their 
broad-band SEDs, in five main classes: spiral, starburst, 
SF-AGN, AGNl and AGN2. 

(ii) computed the rest-frame LPs at 35, 60 and 90 /im 
up to from the 70-, 100- and 160-/xm selected samples 
respectively. 

(iii) integrated the SEDs over Arest =8-1000 /im and com- 
puted the total IR LP up to 2:~4 and studied its evo- 
lution with redshift, finding strong luminosity evolution 
<x{l+zf-^^^'^-^'^ up to 2~1.85, and oc(l-h2)''^"° '^' between 
2~1.85 and 2:~4, combined with a negative density evolu- 
tion oc(l-h2)-°-"±°-'' up to z^l.l and oc(l-h2)-^-^'±''-^* at 
z>l.l and up to z^A. 

(iv) derived the evolution of the comoving total IR lu- 
minosity density, which is found to increase as (l-t-z)^"*" '^ 
up to z~l.l, then to remain nearly constant (decrease as 



(1+z) " ) from 2:~1.1 to 2~2.8, and to decrease as 
(l+z)-«-''±«-^ up to z~4. 

(v) found that the evolution derived for the global IR LP 
is indeed a combination of different evolutionary paths: the 
IR population does not evolve all together "as a whole" , as is 
often assumed in the literature, but is composed of different 
galaxy classes evolving differently: the spiral-SED galaxies 
dominate pm. only at low redshifts {z^ 0.5-0.6), then SF-AGN 



dominate up to z;~2.5, while AGNl and AGN2 start dominating 
the IR luminosity density only at z > 2.5. 

(vi) derived the relative contribution to pm. of MS and 
off-MS sources, which keep nearly constant between z~0.8 
and 2~2.2, with the MS population always dominating. The 
contribution to pia of the off-MS sources shows no significant 
signs of increase with z (from ~19 per cent at 2~0.8-l.25 
to > 22 per cent at 2~1.8-2.2). 

(vii) derived very different evolutionary behaviour in 
terms of different contributions to piR , for galaxies with dif- 
ferent masses, refiecting the downsizing scenario (piR peaks 
at higher redshift with increasing mass). Intermediate-mass 
objects (log(Af/MQ) = 10-ll) always dominate the IR lumi- 
nosity density, increasing with redshift up to 2~1, then re- 
maining nearly constant at higher redshifts (at least up to 
z~2.8), while the contribution of most massive objects in- 
creases even more rapidly with z (at 2~2 it was ~5 times 
higher than today) and continues to grow up to 2~3. 

(viii) described a possible twofold evolutionary scenario 
for IR sources: on the one hand, AGNl and AGN2, representing 
the same population, after an intense starburst phase (due to 
a major merging event, appearing as SF-AGN (SB) SED galax- 
ies), suppress the SF (and shine as AGN) and evolve to red 
massive spheroids; on the other hand, the SF-AGN (Spiral) 
galaxies represent a phase in the life of a star- forming galaxy, 
following a moderate burst of SP (starburst, with SEDs like 
those of local starburst galaxies) and preceding the forma- 
tion of a steady spiral galaxy as we observed in the local 
Universe. 
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Table 3: Rest-frame 35 jim luminosity function 





from the GOODS-S 70-/im sample 
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Table 4: Combined rest-frame 60 jiia luminosity function 
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Table 5: Combined rest-frame 90 yum luminosity function 
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Table 6: PEP total IR luminosity function 



log(LiR/L0) 
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Table 9: PEP total IR lumiiio.sity density 
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1.01±0.71 


1.21±0.91 


1.02±0.52 










log(.1//Mo) 


=10- 


11 


0.99±0.08 


1.65±0.21 


2.14±0.40 


3.54±0.75 


2.94±0.52 


4.05±0.90 


4.18±1.85 


2.28±0.41 


3.18±0.61 


4.47±2.69 


1.48±0.71 


log(Af/Mo) 


=11- 


12 


0.08±0.05 


0.26±0.06 


0.23±0.12 


0.31±0.18 


0.63±0.07 


0.71±0.31 


0.85±0.20 


0.82±0.22 


1.23±0.46 


2.82±2.24 


0.34±0.10 




0.8<z<1.25 


1.25<«<1.8 


1.8<«<2.2 


















on-MS 






5.71±1.01 


4.02±1.27 


7.12±5.03 


















off-MS 






1.41±0.15 


1.19±0.58 


2.01±0.65 



















